Method and apparatus for medical imaging using combined near-infrared optical tomography, fluorescent tomography and ultrasound

ABSTRACT

Methods and apparatus for medical imaging using diffusive optical tomography and fluorescent diffusive optical tomography and ultrasound are disclosed. In one embodiment, the probe comprises emitters and detectors that are inclined at an angle of about 1 to about 30 degrees to a surface of the probe that contacts tissue. In another embodiment, the scanned volume is divided into an inclusion region and a background region. Different voxel sizes are used in the inclusion region and the background region. Appropriate algorithms facilitate a reconstruction of the inclusion region to determine structural and functional features of the inclusion.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims priority to Provisional Application No.60/832,002 filed Jul. 19, 2006, the entire contents of which are herebyincorporated by reference.

STATEMENT OF FEDERALLY FUNDED RESEARCH

This invention was made with support from the United States Governmentunder contract number R01EB002136 awarded by the National Institutes ofHealth and contract number DANK) 17-03-1-0690 awarded by the Departmentof Defense. The Government has certain rights in the invention.

BACKGROUND

This disclosure relates primarily to the field of biological imaging,particularly to medical imaging. More specifically, to medical imagingequipment and methods for medical imaging using combined near-infraredoptical tomography, fluorescent tomography and/or ultrasound.

Diffusive optical tomography (DOT) is a form of computer-generatedtomography wherein near-infrared light (NIR) is directed at a inclusion(e.g., a lesion, a tumor, a cancer, and so forth) and the amount oflight transmitted and/or diffused through the object, and/or reflectedfrom the object, is detected and utilized to reconstruct a digital imageof the target area (e.g., the object can exhibit a differential intransmission and/or diffusion from surrounding tissues). This method ofimaging is of interest for several reasons, for example, differing softtissues exhibit differing absorption, transmission and/or scattering ofnear-infrared light. Therefore, DOT is capable of differentiatingoptical properties between inclusions and normal soft tissues, whereinalternative tomography methods (e.g., Positron Emission Tomography,Magnetic Resonance Imaging, X-Ray, and so forth) cannot. Another exampleis that near-infrared light is non-ionizing to bodily tissues, andtherefore patients can be subjected to repeated light illuminationwithout harm. This in turn allows physicians to increase the frequencyat which they monitor and/or track change in areas of interest (e.g.,lesions, tumors, and so forth). Yet further, due to differences at whichnatural chromophores (e.g., oxygen-hemoglobin) adsorb light differentlyat different wavelengths, optical tomography is capable of supplyingfunctional information such as hemoglobin concentration and oxygensaturation. For these reasons there is much interest in employingoptical tomography for the detection and monitoring of canceroustissues, especially in breast cancer applications.

Although diffusive optical tomography is a promising medical imagingtechnique, DOT imaging methods and DOT apparatus have yet to yield highquality reconstructions of inclusions due to fundamental issues withintense light scattering.

Another method of tomography imaging that is of interest is fluorescentdiffusive optical tomography (FDOT). Fluorescent diffusive opticaltomography is a form of computer-generated tomography wherein anexcitation source (e.g., near-infrared light) is directed at aninclusion labeled by a fluorescence targets or dyes or fluorophores.Upon excitation of the fluorophore, the wavelength of the excitationsource is shifted to a differing wavelength (e.g., a Stokes-shift) as itis emitted by the fluorophore. The emitted light is then detected andutilized to reconstruct a digital image of the target area, which canexhibit a differential in fluorophore concentration from surroundingtissues (e.g., fluorophore take-up). The digital image can be employedto provide functional characteristics about the inclusion, such asvascular endothelial growth factor (VEGF) conjugated with a dyefluorophore. However, FDOT methods have exhibited less than desirablereconstruction accuracy due to imperfect uptake of the fluorophore andbackground fluorophore noise.

While diffusive optical tomography provides benefits over alternativeimaging methods, they suffer from drawbacks. Each of these imagingmethods is confronted with challenges that impede widespread acceptanceand implementation. One of the drawbacks is that the accuracy ofreconstructed hemoglobin concentration and oxygen saturation is lowwithout location information about the inclusion. It is thereforedesirable to develop methods that can provide structural and functionalcharacteristics about inclusions while at the same time facilitatingaccuracy in inclusion location so that accurate information abouthemoglobin concentration and oxygen saturation can be obtained.

SUMMARY

Disclosed herein is a method for medical imaging of an inclusioncomprising imaging a tissue volume with a probe comprising: anultrasound transducer that is operative to provide an on-site estimationof inclusion size and location; a first emitter and a first detector;the first emitter having light of a wavelength of about 400 to about 900nanometers; the first detector detecting light of a wavelength of about400 to about 900 nanometers; a source circuit connected in operationalcommunication to the emitter; a detector circuit connected inoperational communication to the detector; and a central processing unitconnected to the source circuit and the detector circuit; scanning thetissue volume with light having a wavelength of about 400 to about 900nanometers; the tissue volume comprising a first layer and a secondlayer; segmenting the scanned tissue volume into an inclusion regioncomprising a plurality of first voxels and a background regioncomprising a plurality of second voxels; the volume of each second voxelbeing larger than the volume of each first voxel; identifying a tissuelayer thickness and an approximate tilting angle between the tissuelayer and the probe; estimating photon density:φ(r,ω)=f(μ_(a1),μ_(s1)′,μ_(a2),μ_(s2)′), where μ_(a1), μ_(a2), μ_(s1)′,and μ_(s2)′ are absorption and reduced scattering coefficients of firstand second layers, respectively; wherein absorption coefficients havethe subscript “a”, scattering coefficients have the subscript “s”,wherein the subscript “1” represents the first layer, the subscript “2”represents the second layer, and the superscript (′) stands for thereduced scattering coefficient; applying a conjugate gradient method toestimate fitted optical properties of the first layer and the secondlayer; minimizing a least-square objective function given asmin∥φ(r,ω)−φ₀(r,ω)∥²; wherein a forward Jacobian weight matrix

${{Wij} = \left\lbrack {\frac{\partial\varphi_{ij}}{\partial\mu_{aj}},\frac{\partial\varphi_{ij}}{\partial D_{j}}} \right\rbrack},$

that relates the photon density perturbation at detector i and imagingvoxel j with absorption coefficient change Δμ_(aj) and diffusioncoefficient change ΔD_(j), is calculated by using the bulk opticalproperties simulated from the two-layer layer model and given inequation (1) as:

$\begin{matrix}{{\left\lbrack W_{ij} \right\rbrack = \begin{bmatrix}\frac{\partial\varphi_{11}}{\partial\mu_{a\; 1}} & \ldots & \frac{\partial\varphi_{1L}}{\partial\mu_{aL}} & \frac{\partial\varphi_{11}}{\partial D_{1}} & \ldots & \frac{\partial\varphi_{1\; L}}{\partial D_{L}} \\\frac{\partial\varphi_{21}}{\partial\mu_{a\; 1}} & \ldots & \frac{\partial\varphi_{2\; L}}{\partial\mu_{aL}} & \frac{\partial\varphi_{21}}{\partial D_{1}} & \ldots & \frac{\partial\varphi_{2L}}{\partial D_{L}} \\\vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\\frac{\partial\varphi_{M\; 1}}{\partial\mu_{a\; 1}} & \ldots & \frac{\partial\varphi_{ML}}{\partial\mu_{aL}} & \frac{\partial\varphi_{M\; 1}}{\partial D_{1}} & \ldots & \frac{\partial\varphi_{ML}}{\partial D_{L}}\end{bmatrix}},} & (1)\end{matrix}$

where M is the total number of detector readings; L is the total numberof imaging voxels and φ₀(r,ω) is the photon density of the first layer.

Disclosed herein too is a method comprising obtaining an image of atissue volume with a probe comprising an ultrasound transducer that isoperative to provide an on-site estimation of inclusion size andlocation; a first emitter and a first detector; the first emitter havinglight of a wavelength of about 400 to about 900 nanometers; the firstdetector detecting light of a wavelength of about 400 to about 900nanometers; a source circuit connected in operational communication tothe emitter; a detector circuit connected in operational communicationto the detector; and a central processing unit connected to the sourcecircuit and the detector circuit; scanning the tissue volume with lighthaving a wavelength of about 400 to about 900 nanometers; the tissuevolume comprising a first layer and a second layer; segmenting thescanned tissue volume into an inclusion region comprising a plurality offirst voxels and a background region comprising a plurality of secondvoxels; the volume of each second voxel being larger than the volume ofeach first voxel; reconstructing the image using the equation (2)

[U _(sd)]_(M×1) =[W _(L) ,W _(B)]_(M×N) ×[M _(L) ,M _(B)]^(T)_(N×1)  (2) where

W_(L) and W_(B) are weight matrices for the inclusion and backgroundregions, respectively;

[M_(L)] = [∫_(1_(L))Δμ_(a)(r^(′))³r^(′), …  ∫_(N_(L))Δμ_(a)(r^(′))³r^(′)]and[M_(B)] = [∫_(1_(B))Δμ_(a)(r^(′))³r^(′), …  ∫_(N_(B))Δμ_(a)(r^(′))³r^(′)]

are the total absorption distributions of inclusion and backgroundregions, respectively; the weight matrices being calculated based on thebackground absorption coefficient μ _(a) and reduced scatteringcoefficient μ′_(s) measurements obtained from the normal contralateralbreast; a Born approximation being used to relate a scattered fieldU_(sd)(r_(si),r_(di),ω) measured at the source (s) and detector (d) pairi to absorption variations Δμ_(a) (r′) in each volume element of the tworegions within the sample.

Disclosed herein too is an apparatus for medical imaging using diffusiveoptical tomography and fluorescent diffusive optical tomographycomprising; a probe comprising a first emitter and a first detector; thefirst emitter or the first detector being inclined at an angle θ of upto about 30 degrees to a surface of the probe that contacts tissue; theprobe comprising a an ultrasound transducer; a source circuit connectedin operational communication to the emitter; a detector circuitconnected in operational communication to the detector; a centralprocessing unit connected to the source circuit and the detectorcircuit; a display operably connected to the central processing unit;and wherein the central processing unit is capable of processinginformation to provide diffusive optical tomography and fluorescentdiffusive optical tomography.

BRIEF DESCRIPTION OF FIGURES

Referring now to the figures, which are exemplary embodiments, andwherein the like elements are numbered alike.

FIG. 1 is a simplified block diagram of an imaging system;

FIG. 2 is an illustration of the front of a probe, which comprising aplurality of emitters and detectors;

FIG. 3 is an exemplary imaging system;

FIG. 4 (a) is an exemplary illustration of the emission and detection ofradiation using one emitter and two detectors;

FIG. 4( b) is an exemplary depiction of the use of ultrasound and lightfor the detection of shallow inclusions, when the emitters are alignedto be perpendicular to the surface of the tissue;

FIG. 4( c) is an exemplary depiction of the use of ultrasound and lightfor the detection of shallow inclusions, when the emitters are alignedat an angle of about 20 to about 30 degrees to the surface of thetissue;

FIG. 5 is an exemplary schematic diagram of the emitter subsystem;

FIG. 6 is an exemplary schematic diagram of the detector subsystem;

FIG. 7 (a) is a picture illustrating chest-wall induced distortion;

FIG. 7 (b) is an ultrasound image of a small cyst (red arrow) withtitled chest wall marked with an array of white arrows;

FIG. 7 (c) is a plot showing the light reflectance R(ρ) (log(ρ²R(ρ)) vs.ρ, where R(ρ) is measured at the surface of the breast, ρ is the sourcedetector separation. The curve was linear up to about 4 cm ofsource-detector separation marked by the arrow;

FIG. 8( a) shows a B-scan ultrasound image of a normal breast withtilted chest wall marked;

FIG. 8( b) is the corresponding two-layer model generated;

FIG. 9 is an exemplary depiction of the dual-zone mesh based imagereconstruction, the entire tissue volume is segmented based onco-registered ultrasound images into an inclusion region, L (region ofinterest, ROI), and a background region, B.

FIG. 10 is a plot showing the Y % versus the ratio for 4 wavelengthpairs (660 nm/830 nm, 690 nm/830 nm, 730 inn/830 nm and 780 nm/830 nm).

FIG. 11 is a plot showing a reconstructed target optical absorptionversus true values; the circles are results obtained from a two-layermodel and triangles are from a semi-infinite model;

FIG. 12 (a) depicts a B-scan ultrasound image of a two-layer tissuemodel. First layer is made of the homogeneous 0.5% Intralipid solution,and the second layer is made of a silicone tissue model;

FIG. 12 (b) depicts a B-scan ultrasound image of a target located on topof the second layer;

FIG. 12 (c) depicts an absorption map at 780 nm from the semi-infinitemodel;

FIG. 12 (d) depicts a diffusion map at 780 nm from the semi-infinitemodel;

FIG. 12 (e) depicts an absorption map at 780 nm from the layered model;

FIG. 12 (f) depicts a diffusion map at 780 nm from the layered model;

FIG. 13 depicts the co-registered ultrasound images (first column) andreconstructed optical absorption maps of the 3-cm-diameter high-contrastabsorber obtained at 50 MHz (second column), 140 MHz (third column), and350 MHz (fourth column) modulation frequencies;

FIG. 14 depicts the reconstructed optical absorption maps of the3-cm-diameter low-contrast absorber obtained at 50 MHz (first column),140 MHz (second column) and 350 MHz (third column) modulationfrequencies, respectively;

FIG. 15 depicts the co-registered ultrasound images (first column) andreconstructed optical absorption maps of the 1-cm high-contrast absorberobtained at 50 MHz (second column), 140 MHz (third column) and 350 MHz(fourth column) modulation frequencies, respectively;

FIG. 16 depicts the reconstructed optical absorption maps of the 1-cmlow-contrast absorber obtained at 50 MHz (first column), 140 MHz (secondcolumn) and 350 MHz (third column) modulation frequencies; and

FIG. 17 (a) shows the ultrasound image of the cancer, which is locatedbetween 0.5 cm and 3.5 cm in depth;

FIG. 17 (b 1, b 2), (c 1, c 2) and (d 1, d 2) respectively arereconstructed absorption maps obtained from 50 MHz and 140 Mhz at 690nm, 780 nm and 830 nm, respectively of the cancer shown in the FIG. 17(a);

FIG. 17 (e 1, e 2) and (f 1, f 2) respectively depict the totalhemoglobin maps and the blood oxygen saturation maps at 50 MHz and 140MHz;

FIG. 18 (a) is a graphical representation showing amplitude versuseffective reflection coefficient (R_(eff));

FIG. 18 (b) is a graphical representation showing phase versus effectivereflection coefficient (R_(eff)); and

FIG. 19 shows images reconstructed from phantom experimental data: (a)absorbing boundary (R_(eff)≈0), (b) partial reflecting boundary(R_(eff)≈0.6). #1 and #2 correspond to the first slice of x-y image at0.5 cm below the surface and the second slice at 1.0 cm below thesurface, respectively.

DETAILED DESCRIPTION

It will be understood that when an element or layer is referred to asbeing “on,” “interposed,” “disposed,” or “between” another element orlayer, it can be directly on, interposed, disposed, or between the otherelement or layer or intervening elements or layers may be present.

It will be understood that, although the terms first, second, third, andthe like may be used herein to describe various elements, components,regions, layers and/or sections, these elements, components, regions,layers and/or sections should not be limited by these terms. These termsare only used to distinguish one element, component, region, layer orsection from another element, component, region, layer or section. Thus,first element, component, region, layer or section discussed below couldbe termed second element, component, region, layer or section withoutdeparting from the teachings of the present invention.

As used herein, the singular forms “a,” “an” and “the” are intended tocomprise the plural forms as well, unless the context clearly indicatesotherwise. It will be further understood that the terms “comprises”and/or “comprising,” when used in this specification, specify thepresence of stated features, integers, steps, operations, elements,and/or components, but do not preclude the presence or addition of oneor more other features, integers, steps, operations, elements,components, and/or groups thereof.

Unless otherwise defined, all terms (including technical and scientificterms) used herein have the same meaning as commonly understood by oneof ordinary skill in the art to which this invention belongs. It will befurther understood that terms, such as those defined in commonly useddictionaries, should be interpreted as having a meaning that isconsistent with their meaning in the context of the relevant art andwill not be interpreted in an idealized or overly formal sense unlessexpressly so defined herein.

Disclosed herein is a medical imaging apparatus and methods for medicalimaging wherein diffusive optical tomography (DOT) and/or fluorescentdiffusive optical tomography (FDOT) is employed in conjunction withultrasound. In one embodiment, the DOT can be employed in conjunctionwith the FDOT and further in conjunction with ultrasound, xrays,photoacoustic techniques and/or magnetic resonance imaging. In anotherexemplary embodiment, the DOT can be employed in conjunction with theFDOT and further in conjunction with ultrasound. The images can beenhanced using one or more of three algorithms provided below. Theapparatus and methods yield structural information (e.g., position in X,Y, Z coordinates, radius, and so forth) as well as functionalinformation (e.g., fluorophore concentration and hemoglobinconcentration) of an absorbing and fluorescing inclusion within a turbidmedium (e.g., soft tissue).

The method disclosed herein is advantageous in that it can be used forprobing inclusions (e.g., cancers, tumors, lesions, and the like) havinga size greater than 2 to 3 centimeters with a low optical modulationfrequency to ensure deep penetration. The method can also advantageouslybe used to probe smaller inclusions located closer to the skin surfacewith a high optical modulation frequency to ensure high sensitivity.With the real-time co-registered ultrasound guidance on inclusion sizeand location, an optimal light modulation frequency can be selectedwhich may yield more accurate estimates of inclusion angiogenesis andhypoxia.

Diffuse optical tomography (DOT) in the near infrared region (NIR)provides a unique approach for functional based diagnostic imaging.However, the intense light scattering in tissue produced by the DOTdominates the NIR light propagation and makes three-dimensionallocalization of inclusions and accurate quantification of the opticalproperties of these inclusions difficult. Optical tomography guided byco-registered ultrasound (US), magnetic resonance imaging (MRI),photoacoustic techniques or x-ray, has a great potential to overcome theuncertainty of location such inclusions and also to improve lightquantification accuracy.

In this co-registration approach, a region of interest (ROI) can bechosen from a non-optical modality ultrasound or MRI, or photoacoustictechnique, or x-ray to guide DOT on image reconstruction. However, anaccurate delineation of inclusions for DOT from a non-optical modalityis difficult due to different contrast mechanisms. Regions of suspiciousinclusions seen by optical methods could potentially improve theaccuracy of reconstructed optical properties. Described herein is amethod that can optically fine tune the inclusions seen by ultrasound,photoacoustic techniques, MRI or x-rays and then reconstruct thefunctional optical properties of the inclusion, such as opticalabsorption, hemoglobin concentration and fluorescence concentration.

At the outset, it is to be understood that the term inclusion is to beinterpreted as any tissue(s), biological mass(es), biologicalentity(ies), and/or foreign object(s) that can be differentiated fromsurrounding tissue(s), biological mass(es), and/or biologicalentity(ies), using diffusive optical tomography, ultrasound and/orfluorescent diffusive optical tomography. For example, an inclusion canbe a tumor that is disposed within soft tissues, such as a tumor withina female breast, wherein the tumor (e.g., comprising epithelial tissues,masenchymal tissues, and so forth) exhibits dissimilar optical diffusionand fluorophore concentration characteristics from surrounding tissues.Also, the term inclusion as used herein can be used interchangeably withthe terms biological entity and target. Further, the term structuralinformation refers to any information gathered or determined withrespect to the structure, or physical shape, of an inclusion, such asposition (e.g., X, Y, Z coordinates), diameter, mass, volume, shape(e.g., circular, elliptical, and so forth), and so forth. Lastly, theterm functional information is to be interpreted as any informationgathered or determined that can be employed by a physician, operator, orone skilled in the art, to determine additional characteristics aboutthe inclusion.

In one embodiment, frequency domain diffuse optical tomography is usedin combination with an ultrasound scanner for probing inclusions ofvarying size located at different depths under the skin. The frequencydomain diffuse optical tomography uses a plurality of frequencies toprobe inclusions for size and location under the skin. In oneembodiment, the optical tomography uses three frequencies to obtaininformation about

In an exemplary embodiment, the ultrasound scanner is a real-timeco-registered ultrasound scanner that is used to simultaneously providean on-site estimation of inclusion size and location.

In another embodiment, the frequency domain diffuse optical tomographycan be used for accurate characterization of inclusions that that arelocated a depths of less than or equal to about 1 centimeter from theouter surface of the skin. The location of the inclusions is identifiedby real-time ultrasound. This is accomplished by introducing the lightat an angle into the tissue. When light is introduced into a tissue itfollows a banana path, as will be depicted later. By varying the path oflight in the tissue, it is possible to detect inclusions that arelocated at depths of less than or equal to about 1 centimeter from theouter surface of the skin.

Referring now to FIG. 1, an exemplary simplified block diagram of animaging system 200 is illustrated, wherein the imaging system 200comprises a probe 4 that can be disposed on bodily tissue 6 to image aninclusion 8 therein. The probe 4 comprises a first emitter 10 and afirst detector 12, wherein the first emitter 10 is connected inoperational communication to a source circuit 14, and the first detector12 is connected in operational communication to a detector circuit 16.The source circuit 14 and detector circuit 16 are operably connected toa central processing unit 18 (hereinafter referred to as CPU 18), whichis operably connected to a display 20 on which an image of the inclusion8 can be generated. The CPU 18 is capable of controlling the operationof the imaging system 200. The probe also comprises a second emitter 24for emitting and collecting ultrasonic energy into the body tissue 6.

FIG. 2 illustrates a front view of the probe 4 having a plurality offirst emitters 10 and first detectors 12 disposed on a faceplate 30. TheFIG. 2 also depicts a single second emitter (e.g., an ultrasoundtransducer or array) for emitting and collecting ultrasonic energy intothe body tissue. The first emitters 10 are capable of emittingradiation, such as near-ultraviolet light. The first detectors 12 arecapable of detecting radiation emitted by the emitters 10 andfluorescence within the target area (e.g., tissue 6 and inclusion 8).

Any number of first emitters 10 and first detectors 12 can be employedto perforin the imaging function. In an exemplary embodiment, the probe4 can comprise about 1 to about 30 first emitters, specifically about 2to about 20 first emitters and more specifically about 5 to about 10first emitters. A preferred number of first emitters in the probe 4 isabout 9. In another exemplary embodiment, the probe 4 can comprise about1 to about 30 first detectors, specifically about 2 to about 20 firstdetectors and more specifically about 5 to about 10 first detectors. Apreferred number of first detectors in the probe 4 is about 10. It isnoted however that as the number of first emitters 10 and/or firstdetectors 12 increases, the imaging time (e.g., the time elapsed beforeradiation received by a detector can be processed and reconstructed intoan image and displayed on display 20 by CPU 18) can increase due to theadditional information to be processed by the CPU 18.

The first emitters 10 and the first detectors 12 are generally disposedon the surface of the faceplate 30 so that they are in close proximityto the tissue 6 being imaged. The first emitters 10 and first detectors12 can be disposed in any configuration, thereby allowing the imagingvolume to be expanded or localized based on the number and/or spacing offirst emitters 10 and first detectors 12.

In an exemplary embodiment, some of the first emitters 10 can beinclined at an angle θ to the surface of the probe 4, when the probe 4is in its undisturbed form. This is depicted later in the FIG. 4( c).The inclination of the emitters improves image quality for shallowtargets. The angle θ can be about 10 to about 50 degrees, specificallyabout 20 to about 30 degrees as shown in FIG. 4( c).

The probe 4 can also have a plurality of second emitters 24. In anexemplary embodiment, the second emitter is an ultrasound transducer. Inone embodiment, the ultrasound transducer is an ultrasound array. Itwill be recognized that any ultrasound array can be used in the probe 4.For example, the ultrasound array can be 1-dimensional, 2-dimensional,1.5-dimensional or 1.75-dimensional. In an exemplary embodiment, theprobe 4 can have about 1 to about 10 ultrasound transducers. A preferrednumber of ultrasound transducers is 1. In the present embodiment a1-dimensional array is used.

The ultrasound transducer uses ultrasonic energy having a frequencyrange of approximately 20,000 to 10 billion cycles per second (hertz).Ultrasound has different velocities in tissues that differ in densityand elasticity from others. This property permits the use of ultrasoundin outlining the shape of various tissues and organs in the body.

The specific shape of the probe 4 and/or the faceplate 30 is desirablyconfigured to be an ergonomic design that is suited to traverse acrossthe tissue 6 of a patient without causing discomfort to the patient(e.g., the faceplate 30 can comprise rounded edges, a smooth surface,and so forth). The probe can have a completely absorbing surface, i.e.,it can be black in color. However, it can also have a reflectingsurface, i.e., it can be non-black in color. In an exemplary embodiment,it can be white in color.

In addition, the probe 4 can be configured such that it can be hand-heldby an operator. While the exemplary depiction of the probe 4 in the FIG.2 shows a circular cross-sectional area, the cross-sectional area canhave a geometry that is square, rectangular, triangular or polygonal. Inaddition, it is further envisioned the probe 4 can be releasably securedto the conduit (e.g., fiber optic cables, wires, and so forth) thatconnects the probe 4 in operable communication with the source circuit14 and detector circuit 16. Furthermore, the optical fibers deployed onthe probe 4 can be perpendicular to the probe surface or titled withcertain angle to ensure best near infrared light illumination of thebiological tissue under evaluation as shown in FIG. 4( c).

It is desirable for the surface of the probe 4 to be manufactured froman organic polymer, preferably one that is flexible at room temperature,so that it can be used to accommodate the contours of a body whosetissue is under observation. The organic polymer can comprise a widevariety of thermoplastic resins, blend of thermoplastic resins,thermosetting resins, or blends of thermoplastic resins withthermosetting resins. The organic polymer may also be a blend ofpolymers, copolymers, terpolymers, or combinations comprising at leastone of the foregoing organic polymers. The organic polymer can also bean oligomer, a homopolymer, a copolymer, a block copolymer, analternating block copolymer, a random polymer, a random copolymer, arandom block copolymer, a graft copolymer, a star block copolymer, adendrimer, or the like, or a combination comprising at last one of theforegoing organic polymers. Exemplary organic polymers for use in theprobe 4 are elastomers that have glass transition temperatures belowroom temperature.

Examples of the organic polymer are polyacetals, polyolefins,polyacrylics, polycarbonates, polystyrenes, polyesters, polyamides,polyamideimides, polyarylates, polyarylsulfones, polyethersulfones,polyphenylene sulfides, polyvinyl chlorides, polysulfones, polyimides,polyetherimides, polytetrafluoroethylenes, polyetherketones, polyetheretherketones, polyether ketone ketones, polybenzoxazoles,polyphthalides, polyacetals, polyanhydrides, polyvinyl ethers, polyvinylthioethers, polyvinyl alcohols, polyvinyl ketones, polyvinyl halides,polyvinyl nitriles, polyvinyl esters, polysulfonates, polysulfides,polythioesters, polysulfones, polysulfonamides, polyureas,polyphosphazenes, polysilazanes, styrene acrylonitrile,acrylonitrile-butadiene-styrene (ABS), polyethylene terephthalate,polybutylene terephthalate, polyurethane, ethylene propylene dienerubber (EPR), polytetrafluoroethylene, fluorinated ethylene propylene,perfluoroalkoxyethylene, polychlorotrifluoroethylene, polyvinylidenefluoride, or the like, or a combination comprising at least one of theforegoing organic polymers.

Examples of blends of thermoplastic resins includeacrylonitrile-butadiene-styrene/nylon,polycarbonate/acrylonitrile-butadiene-styrene, acrylonitrile butadienestyrene/polyvinyl chloride, polyphenylene ether/polystyrene,polyphenylene ether/nylon, polysulfone/acrylonitrile-butadiene-styrene,polycarbonate/thermoplastic urethane, polycarbonate/polyethyleneterephthalate, polycarbonate/polybutylene terephthalate, thermoplasticelastomer alloys, nylon/elastomers, polyester/elastomers, polyethyleneterephthalate/polybutylene terephthalate, acetal/elastomer,styrene-maleicanhydride/acrylonitrile-butadiene-styrene, polyetheretherketone/polyethersulfone, polyether etherketone/polyetherimidepolyethylene/nylon, polyethylene/polyacetal, or the like.

Examples of thermosetting resins include polyurethane, natural rubber,synthetic rubber, epoxy, phenolic, polyesters, polyamides,polysiloxanes, or the like, or a combination comprising at least one ofthe foregoing thermosetting resins. Blends of thermoset resins as wellas blends of thermoplastic resins with thermosets can be utilized. Anexemplary thermosetting resin is polydimethylsiloxane (PDMS).

Referring now to FIG. 3, an exemplary embodiment of an imaging system200 is illustrated. The imaging system 200 comprises a probe 4 having aplurality of first emitters 10 operably connected to a source circuit 14via optical fibers 32. In one embodiment, the probe has 9 emitters. Thesource circuit 14 comprises an excitation source (ES) 34 that isoptically connected to a primary optical switch (OS1) 36. The excitationsource (ES) 34 comprises multiple excitation elements therein (notshown), such as pigtailed laser diodes capable of emitting near-infraredradiation at 660 nm or 690 nm and near-infrared radiation at 780 nm and830 nm (e.g., commercially available from Thorlabs Inc.) that ismodulated at predetermined frequencies (e.g., 350 MHz, 140 MHz and 50MHz) by an oscillator (OSC2) 90, which is connected thereto.

The primary optical switch (OS1) 36 is capable of selectively connectingthe emissions from any of the excitation elements, or any combination ofexcitation elements, to a secondary optical switch (OS2) 38 (e.g.,commercially available from Piezosystem Jena Inc.). The secondaryoptical switch (OS2) 38 is capable of selectively directing theemissions from the primary optical switch (OS1) 36 connected to anycombination of the nine emitters 10, hence allowing the emission ofradiation through the emitters 10 selected. The primary optical switch(OS1) 36 and the secondary optical switch (OS2) 38 are connected inoperable communication with, and controlled by, CPU 18.

In FIG. 4( a), an exemplary illustration of the emission and detectionof radiation using one first emitter 10 and two first detectors 12. Tobe more specific, first emitter 10 emits radiation 72 through tissue 6upon an inclusion 8 disposed a distance 82 from the first emitter 10.Fluorophores within the inclusion 8 are excited by the radiation 72 andfluoresce in the form of fluorescence signals 74 and 76, which aredetected by a first detector 78 and a second detector 80 that aredisposed a distance 84 and a distance 86 from the first emitter 10,respectively. The fluorophores can be disposed within the inclusion 8via the injection of a fluorophore dye, such as Cy5.5 and the like.

FIG. 4( b) depicts a configuration where the first emitters 10 and thefirst detectors 12 are disposed so as to be parallel to one another. Inthis configuration, the first emitters 10 and the first detectors 12 aregenerally perpendicular to the skin of the patient. In thisconfiguration the light introduced into the tissue travels acrescent-like path, and thus may not impinge on inclusions that arecloser to the skin.

FIG. 4( c) depicts a configuration where the first emitters 10 and/orthe first detectors 12 are inclined at an angle of about 20 to about 30degrees. In this configuration, shallow inclusions located less thanabout 0.5 to about 0.7 centimeters from the surface of the skin can bedetected. Source to detector distances of about 0.5 to about 10centimeters, specifically about 1 to about 8, more specifically about 2to about 7 centimeters enable the detection of shallow inclusions. Aswill be shown later, the probe 4 with a reflection boundary is suitablefor detecting shallow targets with high sensitivity.

FIG. 5 is another exemplary depiction of a source circuit 14. The sourcesubsystem has a plurality of pigtailed laser diodes 34 that emit lighthaving wavelengths of about 400 to about 900 nanometers are in operativecommunication with a corresponding number of radio frequency (RF)oscillators. For each modulation frequency, the outputs of the laserdiodes are sequentially delivered to a plurality of first emitterlocations on the probe through optical switches 36 and 38.

In the embodiment depicted in the FIG. 5, three pigtailed laser diodes34 that emit light having wavelengths of 690 nm, 780 nm and 830 nm (OZOptics Inc) are in operative communication with and three RF oscillatorsof 350 MHz, 140 MHz and 50 MHz. For each modulation frequency, theoutputs of the laser diodes were sequentially delivered to 9 locationson the probe through 3×1 and 1×9 optical switches 36 and 38respectively. The laser diodes were obtained from OZ Optics Inc. and theoptical switches were obtained from PiezoJena Inc.

The ten first detectors 12 on the probe 4 are operably connected to thedetector circuit 16 via optical fibers 52. The detector circuit 16comprises detector sub-circuits 54 for each first detector 12 andoptically connected thereto via portions of optical fibers 52. Eachdetector sub-circuit 54 comprises a collimating system and filter (CSF)54, which is capable of receiving an optical signal (e.g., light) from afirst detector 12, collimating the optical signal, and optionallyfiltering the optical signal to a specific desired frequency range. Theoptical signal emitted from the collimating system and filter (CSF) 54is then directed to a photomultiplier tube (PMT) 56 (e.g., commerciallyavailable from Hamamatsu Inc.) and converted into a voltage, which issubsequently amplified by pre-amp (PA) 58 (e.g., by about 40 mV).

The resulting voltage is mixed with an output carrier signal having apredetermined frequency (e.g., 350.02 MHz, 140.02 MHz and 50.02 MHz) bya local oscillator (OSC1) 60 that is connected in electricalcommunication with the voltage via mixer 62. The heterodyned signalsoutput by mixer 62 are filtered by a narrowband filters (F1) 64 andfurther amplified (e.g., by 30 dB) by amplifier (AMP) 66. The amplifiedsignals are then sampled at a predetermined frequency (e.g., 250 KHz) byan analog to digital conversion (A/D) board inside the CPU 18. Thesignals output by the oscillator (OSC1) 60 are directly mixed with theoutput of oscillator (OSC2) 90 by mixer 68 to produce a reference signal(e.g., a 20 KHz reference signal). The 20 kHz reference signal is thenfiltered by a narrowband filter 70 (e.g., 20 KHz) and provided as inputto the CPU 18.

FIG. 6 is an exemplary depiction of a detector circuit subsystem 54. Thedetection subsystem comprises a plurality of parallel channels. In oneembodiment, the detection subsystem comprises 10 parallel channels. Thelight guides couple the diffusively detected light from the tissue tophoto-multiplier tubes (PMTs) 56. The outputs of the PMTs 56 wereamplified by low-noise preamplifiers 58 and mixed with the correspondingoscillator offset 20 kHz higher to produce both sum and differencesignals. The 20 kHz difference signals were further amplified by 50 dBand bandpass filtered before sampled by A/D converters. Two NationalInstrument (NI) data acquisition cards (8 channels on each card) weresynchronized to acquire data from 10 parallel detection channels. The 20kHz reference signal was produced by directly mixing the outputs of thecorresponding transmit and receive oscillators at each modulationfrequency. The sampled reference signal was used to retrieve phaseinformation. To speed up data acquisition, one or two modulationfrequencies can be selected on site for each breast lesion based on thedepth and target size measured by real-time ultrasound. The entire dataacquisition of 9 10 source-detector positions for three wavelengths wasless than 60 seconds, specifically less than 30 seconds, and morespecifically less than 10 seconds, which was fast enough to acquire datafrom patients and avoid potential motion artifacts.

With reference now once again to the FIG. 3, ultrasound images areavailable from the video output port of any connected commercialultrasound system (US1) 92 and are readily acquired by an imagingcapture card implemented in the CPU. Captured ultrasound images aredirectly sent to the display 20 for review and also for guidingnear-infra-red imaging formation.

The imaging system 200 is capable of diffusive optical tomography (DOT)and ultrasound functions. These imaging techniques are employed toprovide structural information (e.g., size, position, and so forth), aswell as functional information (e.g., fluorophore concentration,hemoglobin concentration, oxygen saturation) in the form of digitalimages of a fluorescing target (e.g., an inclusion) within a turbidmedium that can be viewed on a display 20. The methods for imagingemployed herein are conducted in three main steps, the first is thelocation and size of the inclusion 8 using ultrasound, the second is theestimation of the structural parameters of the inclusion 8, and thethird is the estimation of the functional parameters of the inclusion 8.

As noted above, once an image of an inclusion is obtained using theultrasound transducer and the DOT and/or the FDOT, the followingalgorithms can be used for further determination of the structural andfunctional features of the inclusion. Three sets of algorithms that aredetailed below have been developed, which provide more accurate andreliable optical properties of inclusions of different sizes located atdifferent depths under the skin. As noted above these inclusion sizesand locations are provided by real-time ultrasound images.

The first set of algorithms can be applied to all parts of the body,though it is generally used to segment the imaging volume as breasttissue and chest wall layers and estimate the two layer bulk opticalproperties instead of one layer property to improve backgroundestimation. In short, it can be used to study breast lesions. This firstset of algorithms uses a finite element analysis method to improve theaccuracy of reconstructed lesion optical properties, especially when thechest-wall is located in less than 1.5 to 2 centimeters depths. Theimaging is performed in FEM, which yields improved optical properties oflesions located on top of the chest-wall and closer to the skin surface.

The second set of algorithms reconstructs larger lesions layer by layerin depth instead of the entire volume. This approach significantlyovercomes the ill-conditioned problem in inversion and improves theaccuracy of reconstructed optical properties.

The third set of algorithms described below reconstructs both opticalabsorption and scattering simultaneously in a unique approach that doesnot suffer problems of cross-coupling between estimated absorption andscattering properties.

In another embodiment, the co-registered ultrasound and DOT and/or FDOTcan be used a) to identify the chest-wall depth and the approximatetilting angle between the breast tissue and the chest-wall; b) toestimate the bulk optical properties of both breast tissue andchest-wall layer for accurate weight matrix computation using eitherfinite element method or analytical solution; c) to improve the accuracyof reconstructed lesion optical properties, especially when thechest-wall is located in up to about 5 centimeter depths (where depthsare measured from the under surface of the skin), specifically about 1.0to about 4 centimeter depths, and more specifically about 1.5 to about 2centimeter depths.

In another embodiment, the co-registered ultrasound and DOT and/or FDOTcan be used a) to identify the size of the lesion; b) to segment thelarge lesions of more than 1 centimeters in size into several targetlayers; c) to improve the reconstruction accuracy of lesion opticalproperties by imaging target layer by layer, combined target layers andbackground layers. In yet another embodiment, the co-registeredultrasound and DOT and/or FDOT can be used a) to segment the tissuevolume into target and background regions for simultaneouslyreconstruction of absorption and scattering variations; b) toreconstruct absorption changes using finer and coarse grids andscattering changes using coarse grid only.

Algorithm Set #1:

The first set of algorithms models the breast tissue as a two-layermedium of breast tissue and chest-wall and estimates the opticalproperties of both layers. The first layer depth and angle of theinterface between the two-layers are estimated from the co-registeredultrasound images. The estimated optical properties of both layers areincorporated into the weight matrix calculations. As a result, moreaccurate target optical properties can be reconstructed.

Since conventional ultrasound is used in pulse-echo reflection geometry,it is desirable to acquire optical measurements with the same geometry.Compared with transmission geometry or ring geometry where the lightsources and detectors are deployed on a pair of parallel plates or aring, the reflection geometry has the advantage of probing reducedbreast tissue thickness because patients are scanned in the supineposition. Consequently, inclusions closer to the chest wall can beimaged. In general, the breast tissue thickness has been reduced to lessthan about 3 to about 4 centimeters. However, when the chest wall ispresent within about 1 to about 1.5 centimeters of the skin surface, thesemi-infinite geometry is not a valid assumption for opticalmeasurements and the chest wall underlying the breast tissue affects theoptical measurements obtained from distant source-detector pairs. Ingeneral, the symmetric location of the contralateral normal breast ischosen as the reference site to minimize the effect of the chest wall.

FIG. 7( a) illustrates the problem when multiple source and detectorfibers are deployed on the breast for co-registered imaging. The lightreflectance measurements at distant detector positions from the sourcescomprises mixed signals from both breast tissue and chest wall. Becauseof the superposition of signals from the breast tissue and the chestwall, the signals are distorted. FIG. 7( b) is the ultrasound image ofthe breast where the tilted chest wall marked by white arrows can beseen. FIG. 7( c) plots the light reflectance (log scale) measured at thesurface of the probe versus source-detector separation ρ. If thesemi-infinite assumption were valid, the plot of reflectance (log scale)versus ρ should be linear. However, as the light penetrates deeper intothe tissue and is therefore received by detectors at larger distancesfrom the sources, a deviation from linearity begins to develop which isclearly seen in the plot in FIG. 7( c).

The distorted distant measurements can be quite complex. In somepatients, the measured distant amplitude profiles consist of many randompoints, while in others, the amplitude profiles bend considerably fromexpected linear curves. In general, when the amplitude values are low,the phase profiles behave like random variables. As a result,measurements beyond a certain distance must be removed before imaging.Since the depth of chest wall and average breast tissue absorption andscattering coefficients vary from patient to patient, it is difficult topredict the cutoff source-detector distance. As a consequence, the dataprocessing and image reconstruction must be done off-line for eachindividual patient by examining the data and removing distantmeasurements.

There are many clinical cases in which the measured amplitude and phaseprofiles bend considerably from the expected linear curves due to theinfluence of the shallower chest-wall layer. The resolution betweenthese curved measurements and linearity results a loss of informationobtained during the process. In addition, inclusions may be located ontop of chest walls, and removing curved measurements may partiallyremove the perturbation caused by the inclusion.

In order to resolve these issues and to maximize information obtainedfrom the curved path of the light beams, a simple fitting method basedon a two-layer model can be employed. This model is used to fit themeasured amplitude and phase profiles after filtering out distantsource-detector measurements and to estimate the first and thesecond-layer optical properties, which can be incorporated into theimaging reconstruction. A numerical reconstruction algorithm based onfinite element methods (FEM) has been used for solving complex boundaryconditions. Summaries of these are provided below.

In one embodiment, a 3D cylindrical mesh was generated for forwardcalculation. The radius of the cylinder was large enough to approximatethe semi-infinite geometry, and the height of the cylinder was chosenaccording to the total imaging volume obtained from co-registeredultrasound images. A smooth surface was used to model the tissue and thechest-wall interface. The chest-wall tilting angle α with respect to theprobe and the tissue-chest interface located in depth D was determinedby the co-registered ultrasound images. FIGS. 8( a) and (b) gives anexample of the two-layer model configuration. FIG. 8( a) shows a B-scanultrasound image of a normal breast with tilted chest wall marked andFIG. 8( b) is the corresponding two-layer model generated.

Bulk absorption and reduced scattering coefficients for both first andsecond layers were used to estimate the photon density:φφ(r,ω)=f(μ_(a1),μ_(s1)′,μ_(a2),μ_(s2)′), where μ_(a1), μ_(a2), μ_(s1)′,and μ_(s2)′ are absorption and reduced scattering coefficients of firstand second layers, respectively. Absorption coefficients have thesubscript “a”, while scattering coefficients have the subscript “s”,where the subscript “1” represents the first layer, the subscript “2”represents the second layer, and the superscript (′) stands for thereduced scattering coefficient. To find the best fitting results,reasonably accurate bulk absorption and scattering coefficients μ_(a10)and μ_(s10)′ of the first layer were first obtained by fitting themeasurements from normal breasts after removing distorted distantmeasurements. The initial absorption and scattering coefficients μ_(a20)and μ_(s20)′ of the second layer were set to be the same as those of thefirst layer. Measurements from a homogeneous Intralipid solution, whichhas the same optical properties as the fitted bulk values of the firstlayer, were used as φ₀(r,ω). Intralipid solution is a light scatteringmedium used in photodynamic therapy in both phantom studies and toobtain optimum light distributions from clinical light delivery systems.When the initial estimates of φ₀(r,ω), μ_(a10), μ_(a20), μ_(s10)′, andμ_(s20)′, have been obtained, the conjugate gradient method was appliedto find the best fitted optical properties of the two-layer model byminimizing the least-square objective function given asmin∥φ(r,ω)−φ₀(r,ω)∥².

${{Wij} = \left\lbrack {\frac{\partial\varphi_{ij}}{\partial\mu_{aj}},\frac{\partial\varphi_{ij}}{\partial D_{j}}} \right\rbrack},$

The forward Jacobian weight matrix which relates the photon densityperturbation at detector i and imaging voxel j with absorptioncoefficient change Δμ_(aj) and diffusion coefficient change ΔD_(j), iscalculated by using the bulk optical properties simulated from thetwo-layer layer model and given in equation (1) as:

$\begin{matrix}{{\left\lbrack W_{ij} \right\rbrack = \begin{bmatrix}\frac{\partial\varphi_{11}}{\partial\mu_{a\; 1}} & \ldots & \frac{\partial\varphi_{1L}}{\partial\mu_{aL}} & \frac{\partial\varphi_{11}}{\partial D_{1}} & \ldots & \frac{\partial\varphi_{1\; L}}{\partial D_{L}} \\\frac{\partial\varphi_{21}}{\partial\mu_{a\; 1}} & \ldots & \frac{\partial\varphi_{2\; L}}{\partial\mu_{aL}} & \frac{\partial\varphi_{21}}{\partial D_{1}} & \ldots & \frac{\partial\varphi_{2L}}{\partial D_{L}} \\\vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\\frac{\partial\varphi_{M\; 1}}{\partial\mu_{a\; 1}} & \ldots & \frac{\partial\varphi_{ML}}{\partial\mu_{aL}} & \frac{\partial\varphi_{M\; 1}}{\partial D_{1}} & \ldots & \frac{\partial\varphi_{ML}}{\partial D_{L}}\end{bmatrix}},} & (1)\end{matrix}$

where M is the total number of detector readings and L is the totalnumber of imaging voxels.

Algorithm Set #2

This algorithm provides imaging reconstruction based on a two-layermodel with the layer-structure identified by co-registered ultrasound.

For a large tumor that is larger than two layers in depth (more than 1cm in general), wherein each layer is considered to be about 0.5centimeter, the performance of the dual-zone mesh based inversion isdegraded because of the increased number of fine-mesh grids with unknownoptical properties in the region of interest. As a result, the weightmatrix is more ill-conditioned and underdetermined for inversion. Thereconstruction of large inclusions is improved by implementing thedual-zone mesh one layer (in depth) at a time, combining the targetlayers, and the background layers to obtain the final absorption images.

In the dual-zone mesh based image reconstruction, two types of voxels, afirst voxel and a second voxel, are used. The first voxels lie withinthe boundary of the inclusion, while the second voxels lie outside theboundary of the inclusion (e.g., the background). Those voxels withinthe boundary of the inclusion, are finer than those located outside ofthe boundary of the inclusion. Within the boundary of the inclusion,each first voxel has a volume of about 0.1×0.1×0.1 cubic centimeter(cm³) to about 1×1×1 cm³, specifically about 0.2×0.2×0.2 cm³ to about0.8×0.8×0.8 cm³, to about 0.3×0.3×0.3 cm³ to about 0.6×0.6×0.6 cm³.Within the boundary of the inclusion, each voxel has a preferred volumeof about 0.5×0.5×0.5 cm³.

Outside the boundary of the inclusion, each second voxel has a volume ofabout 0.2×0.2×0.2 cm³ to about 2.0×2.0×2.0 cm³, specifically about1.2×1.2×1.2 cm³ to about 1.8×1.8×1.8 cm³, and more specifically about1.3×1.3×1.3 to about 1.7×1.7×1.7 cm³. Outside the boundary of theinclusion (e.g., for imaging the tissue surrounding the inclusion), eachvoxel has a preferred volume of about 1.5×1.5×1.5 (cm³). The totalimaging volume is chosen to be about 15×15×10 cm³, specifically about12×12×8 cm³, and more specifically about 9×9×4 cm³.

The ratio of voxel sizes inside the boundary of the inclusion (inclusionregion) to that outside the boundary of the inclusion (backgroundregion) is about 3% to about 90%, specifically about 8% to about 50%,more specifically about 10% to about 30%, and even more specificallyabout 12 to about 20%.

With reference now to the FIG. 9, in the dual-zone mesh based imagereconstruction, the entire tissue volume is segmented based onco-registered ultrasound images into an inclusion region, L (region ofinterest, ROI), and a background region, B. A modified Bornapproximation is used to relate the scattered fieldU_(sd)(r_(si),r_(di),ω) measured at the source (s) and detector (d) pairi to absorption variations Δμ_(a)(r′) in each volume element of the tworegions within the sample. The matrix form of the image reconstructionis given by the equation (2)

[U _(sd)]_(M×1) =[W _(L) ,W _(B)]_(M×N) ×[M _(L) ,M _(B)]^(T)_(N×1)  (2)

where W_(L) and W_(B) are weight matrices for the inclusion andbackground regions, respectively.

[M_(L)] = [∫_(1_(L))Δμ_(a)(r^(′))³r^(′), …  ∫_(N_(L))Δμ_(a)(r^(′))³r^(′)]

and

[M_(B)] = [∫_(1_(B))Δμ_(a)(r^(′))³r^(′), …  ∫_(N_(B))Δμ_(a)(r^(′))³r^(′)]

are the total absorption distributions of inclusion and backgroundregions, respectively. The weight matrices are calculated based on thebackground absorption coefficient μ _(a) and reduced scatteringcoefficient μ′_(s) measurements obtained from the normal contralateralbreast.

If large inclusions occupy more than one layer in depth, we divide theregion of interest (ROI) into N zones in N layers denoted as ROI#1,ROI#2, . . . ROI#N respectively. The above dual-zone reconstructionalgorithm will be used one layer at a time and a total of N steps istherefore needed to reconstruct the entire image.

Mathematically, the new reconstruction procedures can be written asshown in the Equations (3), (4) and (5):

Step 1:

[U _(sd)]_(M×1) =[W _(L(ROI#1)) ,W _(B) ₁ ,W _(B) ₁_((other-target-layers)) ,W _(B(BG-layers)) ]×[M _(L(ROI#1)) ,M _(B) ₁,M _(B) ₁ _((other-target-layers)) ,M _(B(BG-layers)) ]T  (3)

Step 2:

[U _(sd)]_(M×1) =[W _(L(ROI#2)) ,W _(B) ₂ ,W _(B) ₂_((other-target-layers)) ,W _(B(BG-layers)) ]×[M _(L(ROI#2)) ,M _(B) ₂,M _(B) ₂ _((other-target-layers)) ,M _(B(BG-layers)) ]T  (4)

. . .

Step N:

[U _(sd)]_(M×1) =[W _(L(ROI#N)) ,W _(B) _(N) ,W _(B) _(N)_((other-target-layers)) ,W _(B(BG-layers)) ]×[M _(L(ROI#2)) ,M _(B)_(N) ,M _(B) _(N) _((other-target-layers)) ,M _(B(BG-layers)) ]T  (5)

where W_(L(ROI#1)), W_(L(ROI#2)), . . . , W_(L(ROI#N)), are weightmatrices of ROI#1, #2 and #N respectively, and W_(B) ₁ , W_(B) ₂ , . . ., W_(B) _(N) are weight matrices of corresponding background regions inthe same depth layer. W_(B) ₁ _((other-target-layers)), W_(B) ₂_((other-target-layers)), W_(B) _(N) _((other-target-layers)), areweight matrices of other target layers excluding the corresponding ROI#target layer, and W_(B(BG-layers)) is the weight matrix of commonbackground layers respectively. After N steps of imaging reconstruction,the total absorption matrix can be obtained as [M_(L(ROI#1)), M_(B) ₁ ,M_(L(ROI#2)), M_(B) ₂ , . . . M_(L(ROI#N)), M_(B) _(N) ,M_(BG(BG-layers))], and the absorption distribution can be computed bydividing the corresponding voxel size in ROIs and in the backgroundregions.

In each reconstruction step, only one layer with finer mesh or ROI # isreconstructed and the other target layers are treated as background.Therefore, the total number of imaging voxels with unknown opticalproperties is always comparable with the total measurements and theinversion is well-conditioned.

Algorithm Set #3

Previously conducted experiments, have shown that when opticalabsorption changes as well as scattering changes are reconstructed, thesimultaneous reconstruction of both absorption and scatteringcoefficients of a target may result in a cross-coupling between thereconstructed absorption and scattering coefficients and this results ina reduced accuracy of the target absorption estimation. The algorithmdescribed below simultaneously reconstructs the target and backgroundabsorption and scattering changes. However, the scattering changes arereconstructed in coarse mesh to significantly reduce the total number ofvoxels with unknown optical properties and to improve the accuracy ofreconstructed absorption coefficients. The absorption estimation isdirectly related to the functional parameter calculation, such ashemoglobin concentration and blood oxygenation.

In the dual-zone mesh based image reconstruction, the entire tissuevolume is segmented based on co-registered ultrasound images into alesion region, L (region of interest, ROI), and a background region, B.A modified Born approximation is used to relate the scattered fieldU_(sd)(r_(si),r_(di),ω) measured at the source (s) and detector (d) pairi to absorption variations Δμ_(a)(r′) in each volume element of tworegions within the sample. The scattered field U_(sd)(r_(si),r_(di),ω)measured at the source (s) and detector (d) pair i is also related tothe diffusion coefficient variations ΔD(r′) in the target andbackground. However, both target and background scattering changes arereconstructed in a coarse grid to significantly reduce the number ofunknowns and improve the accuracy of reconstructed optical absorptionchanges. The matrix form of image reconstruction is given by theEquation (6)

[U _(sd)]_(M×1) =[W _(L) ,W _(B) ,W _(S)]_(M×N) ×[M _(L) ,M _(B) ,M_(S)]^(T) _(N×1)  (6)

where W_(L) and W_(B) are absorption weight matrices for lesion andbackground regions, respectively. W_(S) is the scattering weight matrixfor the entire imaging volume.

[M_(L)] = [∫_(1_(L))Δμ_(a)(r^(′))³r^(′), …  ∫_(N_(L))Δμ_(a)(r^(′))³r^(′)]and[M_(B)] = [∫_(1_(B))Δμ_(a)(r^(′))³r^(′), …  ∫_(N_(B))Δμ_(a)(r^(′))³r^(′)]

are the total absorption distributions of inclusion and backgroundregions, respectively.

[M_(S)] = [∫₁Δ D(r^(′))³r^(′), …  ∫_(N)Δ D(r^(′))³r^(′)]  

is the total diffusion distribution of the entire tissue volumesegmented with the coarse grid. The weight matrices are calculated basedon the background absorption coefficient μ _(a) and reduced scatteringcoefficient μ′_(s) measurements obtained from the normal contralateralbreast.

By assuming that the major chromophores are deoxygenated (deoxyHb) andoxygenated (oxyHb) hemoglobin in the wavelength range studied, a directestimate of the total hemoglobin concentrationtotalHb(r′)=deoxyHb(r′)+oxyHb(r′) and oxygenation saturation

${Y\; \%} = {\frac{{{oxy}{Hb}}\left( r^{\prime} \right)}{{{{oxy}{Hb}}\left( r^{\prime} \right)} + {{deoxy}\left( r^{\prime} \right)}}100\%}$

can be obtained from the Equation (7) and (8):

$\begin{matrix}{{{{total}{Hb}}\left( r^{\prime} \right)} = {{\frac{1}{\Delta}\left\{ {{\left\{ {ɛ_{{HbO}_{2}}^{\lambda_{2}} - ɛ_{Hb}^{\lambda_{2}}} \right\} {\mu_{a}^{\lambda_{1}}\left( r^{\prime} \right)}} + {\left\{ {ɛ_{Hb}^{\lambda_{1}} - ɛ_{{HbO}_{2}}^{\lambda_{1}}} \right\} {\mu_{a}^{\lambda_{2}}\left( r^{\prime} \right)}}} \right\}}\mspace{146mu} = {{{{weight}\left( \lambda_{1} \right)}{\mu_{a}^{\lambda_{1}}\left( r^{\prime} \right)}} + {{{weight}\left( \lambda_{2} \right)}{\mu_{a}^{\lambda_{2}}\left( r^{\prime} \right)}}}}} & \begin{matrix}(7) \\(8)\end{matrix} \\{{{{and}\mspace{14mu} Y\%} = {\frac{{{- ɛ_{Hb}^{\lambda_{2}}}\frac{\mu_{a}^{\lambda_{1}}\left( r^{\prime} \right)}{\mu_{a}^{\lambda_{2}}\left( r^{\prime} \right)}} + ɛ_{Hb}^{\lambda_{1}}}{{\left( {ɛ_{{HbO}_{2}}^{\lambda_{2}} - ɛ_{Hb}^{\lambda_{2}}} \right)\frac{\mu_{a}^{\lambda_{1}}\left( r^{\prime} \right)}{\mu_{a}^{\lambda_{2}}\left( r^{\prime} \right)}} - \left( {ɛ_{{HbO}_{2}}^{\lambda_{1}} - ɛ_{Hb}^{\lambda_{1}}} \right)}100\%}},} & (9)\end{matrix}$

where μ_(a) ^(λ) ¹ (r′) and μ_(a) ^(λ) ² (r′) are optical absorption atvolume element r′ and Δ=ε_(Hb) ^(λ) ¹ ε_(HbO) ₂ ^(λ) ² −ε_(HbO) ₂ ^(λ) ¹ε_(Hb) ^(λ) ² .

In principle, any two wavelengths in the NM window can be used tocompute the total hemoglobin concentration from equation (7). However,positive weights of μ_(a) ^(λ) ¹ (r′) and μ_(a) ^(λ) ² (r′) in equation(7) ensure that the estimated total hemoglobin concentration at eachvoxel is positive. This is useful in the background region where theabsorption coefficients are small and consequently the relativeestimation errors of absorption coefficients are larger than those inthe lesion region.

Similarly, any two wavelengths in the NIR window can be used to computethe oxygenation saturation from equation (4). However, the sensitivityof Y % to ratio

$\frac{\mu_{a}^{\lambda_{1}}\left( r^{\prime} \right)}{\mu_{a}^{\lambda_{2}}\left( r^{\prime} \right)}$

is significantly different for different wavelength pairs. FIG. 10 is aplot showing the Y % versus the ratio for four wavelength pairs (660/830nm, 690/830 nm, 730/830 nm and 780/830 nm). For the wavelength pair660/830 nm, the Y % lies in the physiological range between 0.3 and 4.4

$\left( {{e.g.},\mspace{14mu} {0.3 < \frac{\mu_{a}^{\lambda_{1}}\left( r^{\prime} \right)}{\mu_{a}^{\lambda_{2}}\left( r^{\prime} \right)} < 4.4}} \right).$

However, for the wavelength pair of (780 nm and 830 nm), the Y % isgreater than 100% for

$\frac{\mu_{a}^{\lambda_{1}}\left( r^{\prime} \right)}{\mu_{a}^{\lambda_{2}}\left( r^{\prime} \right)}$

is less than or equal to about 0.7 and reduces quickly to zero for aratio greater than or equal to about 1.4. Any estimation error on μ_(a)^(λ) ¹ (r′) and μ_(a) ^(λ) ² (r′) can cause the estimated Y % to falloutside the expected physiological range, particularly when the (r′) andμ_(a) ^(λ) ² (r′) are small in the background tissue regions. Forexample, the estimated absorption coefficients from normal breasts arewithin 0.01 to 0.06 cm⁻¹ at 780 nm and 0.02 to 0.07 cm⁻¹ at 830 nm. Thecomputed blood oxygen saturation values using these backgroundabsorption coefficients may exceed 100%. Considering both factors,robust estimation of total hemoglobin concentration favors thewavelength pair of 780 nm and 830 nm whereas blood oxygenationsaturation estimation is best for the wavelength pair of 660 nm and 830nm.

With the real-time co-registered ultrasound guidance on inclusion sizeand location, an optimal light modulation frequency can be selectedwhich may yield more accurate estimates of lesion angiogenesis andhypoxia. As will be shown in the examples below, large and smallabsorbers of both high and low optical contrasts have shown that a highmodulation frequency, such as 350 MHz, is preferable for probing smalllesions closer to the surface while a low modulation frequency, such as50 MHz, is desirable for imaging deeper and larger lesions.

Data from many clinical cases have shown that NIR measurements at 660 nmare subject to higher background light scattering because the reducedscattering coefficient μ_(s)′=aλ^(−b) increases with the decreasedwavelength. In this prototype system, we have chosen the wavelength pairof (690 nm, 830 nm) for blood oxygen estimation. 690 nm is the longestwavelength in the visible range and laser diodes at this wavelength areavailable. As a result, three typical wavelengths of 690 nm, 780 nm, and830 nm are implemented in our system for estimation of total hemoglobinconcentration and blood oxygen saturation.

The present invention can be embodied in the form ofcomputer-implemented processes and apparatuses for practicing thoseprocesses. The present invention can also be embodied in the form ofcomputer program code containing instructions embodied in tangiblemedia, such as floppy diskettes, CD-ROMs, hard drives, or any othercomputer-readable storage medium, wherein, when the computer programcode is loaded into and executed by a computer, the computer becomes anapparatus for practicing the invention. The present invention can alsobe embodied in the form of computer program code, for example, whetherstored in a storage medium, loaded into and/or executed by a computer,or transmitted over some transmission medium, such as over electricalwiring or cabling, through fiber optics, or via electromagneticradiation, wherein, when the computer program code is loaded into andexecuted by a computer, the computer becomes an apparatus for practicingthe invention. When implemented on a general-purpose microprocessor, thecomputer program code segments configure the microprocessor to createspecific logic circuits.

The invention is further illustrated by the following non-limitingexamples.

EXAMPLES Example 1

This example was conducted to validate the accuracy of the fittingmethod for estimating the two-layer optical properties. Towell-condition the inversion, a dual-zone mesh method with a targetregion (region of interest (ROI)) identified by co-registered ultrasoundand a background region for weight matrix calculation was used. Forsimulation and phantom experiments, the imaging voxel size within ROIwas chosen as 0.25 cm×0.25 cm×0.5 cm, the pixel size outside of ROI waschosen as 0.5 cm×0.5 cm×1 cm. For clinical cases, the inclusion boundarywas not well identified in general and the ROI was chosen to be muchlarger than the inclusion boundary measured by ultrasound. Because ofthe dual-mesh scheme, the total number of unknown voxels wassignificantly reduced.

Imaging experiments were performed with simulated targets imbedded in atwo-layer media and reconstructed with the correspondinglayered-structure using a semi-infinite model. In these imagingexperiments performed with simulated targets, the chest wall layers wereplaced at depth of 1.3 cm, 1.6 cm and 2.0 cm, respectively. The firstlayer absorption and scattering parameters (μ_(a1) and μ_(s1)′) werekept as 0.02 cm⁻¹ and 6 cm⁻¹ respectively, and the second layerabsorption and scattering parameters (μ_(a2) and μ_(s2)′) were kept at0.08 cm⁻¹ and 10 cm⁻¹, respectively.

For each set of simulations, the optical absorption coefficient of thetarget placed at [0, 0, 1.0 cm] was varied in an amount of 0.1 cm⁻¹ to0.25 cm⁻¹. When the chest wall was located at 1.6 cm and 2.0 cm, thetarget was a 1 cm cube. When the chest wall was located at 1.3 cm, asmaller target of size 1 cm×1 cm×0.5 cm was used. The reduced scatteringcoefficient of the target was always kept at 15 cm⁻¹. The reconstructedaverage μ_(a) obtained from the corresponding two-layer model and thesemi-infinite model are shown in the FIG. 11. The average was computedwithin 6 dB of the maximum value. Results obtained from the layeredmodel are more accurate than those obtained from the semi-infinitemodel. The improvement of reconstructed μ_(a) varies from 10.6% to 13.8%when the second layer is located at 1.3 cm depth. As the first layerthickness increases to 2 cm, the improvement is negligible.

To evaluate the simulation results, two sets of experiments wereconducted. The top layer was made of homogeneous 0.5% Intralipidsolution of fitted values μ_(a)=0.02 cm⁻¹ and μ_(s)′=4.78 cm⁻¹, and thebottom layer was a silicone having values μ_(a)=0.08 cm⁻¹ and μ_(s)′=8.3cm⁻¹. The silicone was placed 1.4 cm and 1.8 cm beneath the Intralipidsolution in the first and second experiment set, respectively.Measurements were made from the homogeneous Intralipid solution, thelayered medium, and the layered medium with a 1 cm³ silicone absorber.Two absorbers of high and low optical contrasts were used. Thecalibrated μ_(a) of the high contrast absorber was between 0.2 cm⁻¹ and0.3 cm⁻¹, and the μ_(s)′ value was 9.6 cm⁻¹. The calibrated μ_(a) of thelow contrast absorber was 0.07 cm⁻¹, and the μ_(s)′ was 9.6 cm⁻¹.

FIG. 12 (a) shows a B-scan (background scan) ultrasound image of atwo-layer medium for reference measurement. The center depth of thebottom layer is 1.4 cm underneath the Intralipid surface. FIG. 12 (b) isa B-scan ultrasound image of the target located on top of the siliconephantom. The high absorber is located at [0, 0, 0.9 cm] inside theIntralipid solution. With the reported fitting procedure, the estimatedbulk absorption and reduced scattering coefficients for the top layerwere 0.025 cm⁻¹ and 6.6 cm⁻¹, and the corresponding values for thesecond layer were 0.078 cm⁻¹ and 7.6 cm⁻¹. FIG. 12 (c) and (d) are thereconstructed absorption and diffusion maps at 780 nm obtained from thesemi-infinite model. FIG. 12 (e) and (f) are the reconstructedabsorption and diffusion maps obtained from the layered model at thesame wavelength. The uniformity of the images is slightly improved byapplying the layered model. Table I gives a comparison of thereconstructed results obtained from the semi-infinite model and thelayered model. When the semi-infinite model was used for reconstruction,the average reconstructed absorption coefficients within the top andbottom layers of the target region were 0.147 cm⁻¹ and 0.147 cm⁻¹,respectively, while the average reconstructed μ_(s)′s within the top andbottom layers of the target region were 11.6 cm⁻¹ and 13.1 cm⁻¹,respectively. When the layered model was used, the average reconstructedabsorption coefficients of the top and bottom layers of the targetregion were 0.17 cm⁻¹ and 0.16 cm⁻¹, respectively, while the averagereconstructed μ_(s)′s of the top and bottom layers of the target regionwere 12.2 cm⁻¹ and 11.4 cm⁻¹, respectively. The accuracy ofreconstructed average μ_(a) and μ_(s)′ was improved by 9.4% and 5.6%assuming μ_(a)=0.2 cm⁻¹, respectively, when the corresponding two-layermodel instead of the semi-infinite model was used.

The reconstruction results obtained from the low contrast absorbercentered at [0, 0, 0.9 cm] with the second silicone layer located 1.4 cmdeep in the Intralipid solution are given in Table I. Using the reportedfitting procedure, the estimated bulk absorption and reduced scatteringcoefficients of the top layer were 0.025 cm⁻¹ and 6.6 cm⁻¹, and thecorresponding values of the second layer were 0.082 cm⁻¹ and 8.3 cm⁻¹.When the semi-infinite model was used for reconstruction, the averagereconstructed absorption coefficients within the top and bottom layersof the target region were 0.068 cm⁻¹ and 0.065 cm⁻¹, respectively, whilethe average reconstructed μ_(s)′s within the top and bottom layers ofthe target region were 8.56 cm⁻¹ and 8.6 cm⁻¹, respectively. When thelayered model was used for reconstruction, the average reconstructedabsorption coefficients within the top and bottom layers of the targetregion were 0.074 cm⁻¹ and 0.078 cm⁻¹, respectively, while the averagereconstructed μ_(s)′s within the top and bottom layers of the targetregion were 8.54 cm⁻¹ and 8.96 cm⁻¹, respectively. The reconstructedaverage values of μ_(a) and μ_(s)′ were very similar when thecorresponding two-layer model instead of the semi-infinite model wasused.

The results of second set experiment are also given in Table I. When thesecond layer was located at 1.8 cm deep in the intralipid, theimprovement of using the layered model compared with using semi-infinitemodel was smaller than that obtained from the first set of experiments.For the high contrast absorber, the reconstructed absorptioncoefficients within the top and bottom layers of the target region wereimproved from 0.187 cm⁻¹ and 0.173 cm⁻¹ to 0.198 cm⁻¹ and 0.188 cm⁻¹,respectively, while the average reconstructed μ_(s)′s within the top andbottom layers of the target region were changed from 9.1 cm⁻¹ and 9.09cm⁻¹, to 9.02 cm⁻¹ and 8.99 cm⁻¹, respectively.

By applying the two-layer model instead of the semi-infinite model, theaccuracy of reconstructed top layer μ_(a) and bottom layer μ_(s)′ wasimproved by 6.5% of the expected value assuming μ_(a)=0.2 cm⁻¹, whileμ_(s)′ changed slightly. For the low contrast absorber, slight changewas observed between the results obtained from the layered model and thesemi-infinite model. When the semi-infinite model was used, thereconstructed absorption coefficients within the top and bottom layersof the target region were 0.073 cm⁻¹ and 0.076 cm⁻¹, respectively, whilethe average reconstructed μ_(s)′ within the top and bottom layers of thetarget region were 10.26 cm⁻¹ and 9.16 cm⁻¹, respectively. When thelayered model was used, the reconstructed absorption coefficients withinthe top and bottom layers of the target region were 0.074 cm⁻¹ and 0.075cm⁻¹, respectively, while the average reconstructed μ_(s)′s within thetop and bottom layers of the target region were 10.24 cm⁻¹ and 9.07cm⁻¹, respectively.

TABLE 1 High contrast absorber Low contrast absorber Target located ontop of the second layer phantom μ_(a) μ′_(s) μ_(a) μ′_(s) Expectedvalues (cm⁻¹) 0.2-0.3 9.6 0.07 9.6 Second layer Reconstructed0.147/0.147 11.60/13.13 0.068/0.065  8.56/8.60 at 1.4 cm values* using asemi-infinite model Reconstructed 0.172/0.160 12.20/11.39 0.074/0.078 8.54/8.96 values* using a layered model Second layer Reconstructed0.187/0.173 9.10/9.09 0.073/0.076 10.26/9.16 at 1.8 cm values* using asemi-infinite model Reconstructed 0.198/0.188 9.02/8.99 0.074/0075 10.24/9.07 values* using a layered model *Mean values obtained withintarget region for each target layer.

Thus, with the proposed fitting method, optical properties of thetwo-layer model can be obtained with acceptable accuracy. The accuracyof μ_(s)′ for the second layer can always reach about 105% to about114.6% of the estimated value, and the accuracy of corresponding μ_(a)is around about 82.7% to about 160% of the estimated value. The accuracyof μ_(a) for the first layer is about 70% of the estimated value, andthe accuracy of corresponding μ_(s)′ is around about 130% to about 137%of the estimated value. As expected, the greater thickness of the firstlayer results in greater errors in the derived optical coefficients.

Example 2

This example was conducted to evaluate the possible effect of fine-meshzone selection on accuracy of the reconstructed absorption coefficients.In order to do so, the target was reconstructed using regions ofinterest (ROIs) twice the size of the real target in x-y dimensions andthree times the size of the target in x-y dimensions. The results showthat the ROI has negligible effect on reconstructed optical propertiesof larger 3-centimeters absorbers. For smaller 1-centimeters absorbers,it is desirable for the region of interest to be twice the size of theabsorber to ensure accurate reconstruction of target optical properties.

Absorbers having diameters of 1-cm and 3-cm of low contrast (calibratedvalue at 780 nm μ_(a)=0.07 cmP^(−1P), μ_(s)′=5.50 cmP^(−1P)) and highcontrast (calibrated value at 780 nm μ_(a)=0.23 cmP^(−1P), μ_(s)′=5.45cmP^(−1P)) were imbedded in Intralipid™ solution of 0.07% concentration.The fitted background μ_(a) and μ_(s)′ were in the range of 0.02 to0.031 cmP^(−1P) and 4.9 to 6.1 cmP^(−1P) in different sets ofexperiments. In each set of 1-cm absorber experiments, the absorber wascentered at depths of 0.5, 1.0, 1.5, 2.0, 2.5, 3, and 3.5 cm from thecenter of the probe. In each set of 3-cm absorber experiments, theabsorber was centered at depths of 2, 2.5, 3, and 3.5 cm. For eachtarget location, three sets of measurements obtained from threemodulation frequencies were used to reconstruct target absorption maps.The average μ_(a) value measured within full width and half maximum(FWHM) of each x-y target layer was used to quantitatively compare thereconstruction results.

To evaluate the possible effect of region of interest or fine-mesh zoneselection on the accuracy of the reconstructed absorption coefficients,the target using a region of interest twice the size of the real targetin x-y dimensions (2×) and three times (3×) the size of the target inx-y dimensions was reconstructed. Results obtained from these two ROIswere compared.

Clinical studies were performed at the Radiology Department of theUniversity of Connecticut Health Center. The Health Center IRB committeeapproved the human subject protocol. A state-of-the art US machine, ATLSono CT, was used for the study. Ultrasound images and opticalmeasurements were acquired simultaneously at multiple locationsincluding the lesion region, a normal region of the same breast, and asymmetric normal region of the contralateral breast. The optical dataacquired at normal regions were used as reference to calculate thescattered field caused by inclusion.

FIG. 13 shows the co-registered ultrasound images (first column) andreconstructed optical absorption maps of the 3-cm-diameter high-contrastabsorber obtained at 50 MHz (second column), 140 MHz (third column), and350 MHz (fourth column) modulation frequencies. The reconstructedabsorption maps shown in FIG. 13 were obtained at 780 nm with the ROI inx-y dimensions equal to twice the size of the real target. The sametarget was translated in 0.5 cm depth steps to target centers at 2 cm,2.5 cm, 3 cm, and 3.5 cm. In other words, slices at distances of 0.5centimeter were examined for the image. In each image, there are 9 x-yimages of 9 cm by 9 cm obtained from 0.5 cm (slice #1) to 4.5 cm (slice#9) in depth with steps of 0.5 cm. The slices are numbered from left toright and top to bottom.

Ideally, the target should uniformly occupy 6 layers in z direction.However, the last two target layers close to the bottom were visibleonly in images obtained at 50 Mhz when the target center was located at2.0, 2.5 and 3.0 cm, respectively. In images obtained at 140 Mhz, thetop four target layers were seen whereas in images obtained at 350 Mhzonly the top three target layers were seen. This is due to thesignificant light absorption of the top target layers that limits thenumber of photons penetrating into the deep layers and ultimatelydetected at the surface. This light-shadowing phenomenon is very similarto the well-known ultrasound shadow effect of a large inclusion, wherebottom layers of the inclusion or normal tissues underneath the lesionappear black in the image (little sound reflection) due to significantsound absorption of top layers. The shadowing effect is more pronouncedat high RF light modulation frequencies.

Table 2 provides the reconstructed average target μ_(a) measured withinFWHM of each target layer and obtained at three modulation frequencies.The first column lists the target center depth and the second throughfourth columns present the reconstructed average μ_(a) obtained at thecorresponding layer for 50, 140, and 350 MHz, respectively. Results arelisted for ROIs of two (2×) and three (3×) the real target size in x-ydimensions. The percentage ratio of the calculated to the calibratedvalue of μ_(a)=0.23 cmP^(−1P) is given in parenthesis. The average μ_(a)calculated from the first five target layers is also given at the bottomof each row. The highest reconstructed values, (68%, 50 Mhz) and (63%,140 Mhz), are obtained at the third target layer when the target centeris located at 2 cm depth. The average values over the first five targetlayers are 59% (50 Mhz) and 47% (140 Mhz), respectively. If the data areaveraged over the first four layers, the reconstruction accuracy at 140Mhz improves to 53%. For 350 Mhz, the highest value, 83%, is obtained atthe second target layer and the average value is 48%. If the data areaveraged over the first three layers, the reconstruction accuracyimproves to 67%. The reconstructed average values drop quickly when thetarget center is translated deeper. In our clinical study, the breast isscanned with the combined probe compressed against the chest wall. Thelesions are often located within a 3 to 4 cm depth. Clinical imagingconditions therefore most closely resemble targets at the 2 cm or 2.5 cmdepths.

TABLE 2 Center Depth 50 MHz 140 MHz 350 MHz High Contrast ROI ROI ROI: 3cm absorber Layers 2 x 3 x Layers 2 x 3 x Layers 2 x 3 x 2.0 cm #1:0.123(54) 0.118(51) #1: 0.110(48) 0.107(47) #1: 0.140(61) 0.139(60) #2:0.135(59) 0.132(57) #2: 0.123(54) 0.121(53) #2: 0.192(83) 0.190(83) #3:0.156(68) 0.153(67) #3: 0.144(63) 0.143(62) #3: 0.149(65) 0.148(64) #4:0.152(66) 0.150(65) #4: 0.111(48) 0.111(48) #4: 0.040(18) 0.036(16) #5:0.110(48) 0.108(47) #5: 0.055(24) 0.055(24) #5: 0.029(13) 0.029(13) #6:0.073(31) 0.071(31) #6: 0.024(11) 0.023(10) #6: 0.027(12) 0.027(12) Ave:0.135(59) 0.132(57) Ave 0.109(47) 0.107(47) Ave: 0.110(48) 0.108(47) 2.5cm #1: 0.123(53) 0.120(52) #1: 0.109(47) 0.107(47) #1: 0.159(69)0.159(69) #2: 0.143(62) 0.140(61) #2: 0.132(57) 0.130(57) #2: 0.178(77)0.177(77) #3: 0.143(62) 0.139(60) #3: 0.112(48) 0.110(48) #3: 0.123(53)0.122(53) #4: 0.110(48) 0.109(47) #4: 0.067(29) 0.066(29) #4: 0.043(19)0.038(17) #5: 0.074(32) 0.072(31) #5: 0.038(17) 0.038(17) #5: 0.029(13)0.028(12) #6: 0.048(21) 0.044(19) #6: 0.023(10) 0.023(10) #6: 0.026(11)0.026(11) Ave: 0.119(52) 0.116(50) Ave: 0.092(40) 0.090(39) Ave:0.106(46) 0.105(46) 3.0 cm #1: 0.119(52) 0.117(51) #1: 0.107(47)0.106(46) #1: 0.165(72) 0.163(71) #2: 0.126(55) 0.123(53) #2: 0.109(47)0.107(47) #2: 0.167(73) 0.167(73) #3: 0.103(45) 0.102(44) #3: 0.083(36)0.081(35) #3: 0.088(38) 0.088(38) #4: 0.075(33) 0.073(32) #4: 0.057(25)0.057(25) #4: 0.034(15) 0.030(13) #5: 0.050(22) 0.047(20) #5: 0.032(14)0.028(12) #5: 0.026(11) 0.026(11) #6: 0.040(18) 0.036(16) #6: 0.024(11)0.023(10) #6: 0.025(11) 0.025(11) Ave: 0.095(41) 0.092(40) Ave:0.078(34) 0.076(33) Ave: 0.096(42) 0.095(41) 3.5 cm #1: 0.095(41)0.092(40) #1: 0.090(39) 0.087(38) #1: 0.141(61) 0.140(61) #2: 0.086(37)0.083(36) #2: 0.091(39) 0.088(38) #2: 0.110(48) 0.110(48) #3: 0.066(29)0.064(28) #3: 0.076(33) 0.075(33) #3: 0.042(18) 0.042(18) #4: 0.049(21)0.043(19) #4: 0.048(21) 0.048(21) #4: 0.026(11) 0.026(11) #5: 0.040(17)0.037(16) #5: 0.030(13) 0.027(12) #5: 0.025(11) 0.025(11) Ave: 0.067(29)0.064(28) Ave: 0.067(29) 0.065(28) Ave: 0.069(30) 0.069(30)

This controlled study clearly shows that large absorbing tumors areunder-reconstructed and advanced system and non-linear algorithmicdevelopments are needed to improve the reconstruction accuracy of targetoptical properties. CW systems offer the best penetration but will notprovide phase information that enhances sensitivity for imaging smallerlesions.

One issue is the selection of ROI or fine-mesh zone used in imagereconstruction. As can be seen from Table I, the ROI has negligibleeffect as the maximum change observed from two different ROIs (2× and3×) is only 3% for targets located at different depths and probed withdifferent modulation frequencies.

FIG. 14 shows the reconstructed optical absorption maps of the3-cm-diameter low-contrast absorber obtained at 50 Mhz (first column),140 Mhz (second column) and 350 Mhz (third column) modulationfrequencies, respectively. Since co-registered ultrasound images of thelow-contrast absorber are essentially the same as the images ofhigh-contrast absorber shown in FIG. 13, these images are not shown inFIG. 14. Table 3 provides the reconstructed average target μ_(a)obtained at three modulation frequencies. The 3× data was similar to the2× data and therefore only the averages are shown in the Table 3. Thepercentage ratio of the calculated value to the calibrated value ofμ_(a)=0.07 cmP^(−1P) is given in parenthesis. It is clear that theaccuracy of reconstructed target absorption maps has been improved inall modulation frequencies. At 50 Mhz, the average reconstructed valueover the first five-layers is 136%, 124%, 105%, and 87% for targetslocated at 2 cm, 2.5 cm, 3 cm and 3.5 cm, respectively. At 140 Mhz, theaverage reconstructed value over the first five-layers is 101%, 90%,81%, and 69%, for the corresponding target locations. If the data areaveraged over the first four layers, the reconstruction accuracyimproves to 113%, 102%, 90%, and 76%, respectively. At 350 Mhz, thetarget mass only appears at top two to three layers. The averagereconstructed value over the first five-layers is 86%, 64%, 60% and 29%,at the corresponding target locations. The average reconstructed valueconsidering the first three-layers is 104%, 86%, 90%, and 41% for therespective depths.

TABLE 3 Center-Depth Low-Contrast 50 Mhz 140 Mhz 350 Mhz 3-cm absorberROI 2 x 3 x ROI: 2 x 3 x ROI: 2 x 3 x   2 cm #1: 0.090(129) #1:0.074(106) #1: 0.073(104) #2: 0.097(139) #2: 0.082(117) #2: 0.103(146)#3: 0.106(151) #3: 0.093(133) #3: 0.044(63) #4: 0.103(146) #4: 0.067(96)#4: 0.049(69) #5: 0.080(114) #5: 0.039(55) #5: 0.033(48) #6: 0.057(81)#6: 0.025(36) #6: 0.017(25) Ave: 0.095(136) 0.089(127) Ave: 0.071(101)0.071(101) Ave: 0.06(86) 0.059(85) 2.5 cm #1: 0.090(129) #1: 0.077(110)#1: 0.076(108) #2: 0.100(143) #2: 0.089(127) #2: 0.049(70) #3:0.099(142) #3: 0.072(103) #3: 0.056(80) #4: 0.083(118) #4: 0.046(66) #4:0.030(43) #5: 0.061(87) #5: 0.032(45) #5: 0.013(18) #6: 0.046(66) #6:0.025(36) #6: 0.007(10) Ave: 0.087(124) 0.083(119) Ave: 0.063(90)0.060(86) Ave: 0.045(64) 0.044(63) 3.0 cm #1: 0.087(125) #1: 0.074(106)#1: 0.095(136) #2: 0.091(130) #2: 0.073(104) #2: 0.066(94) #3:0.080(114) #3: 0.060(86) #3: 0.027(39) #4: 0.062(89) #4: 0.044(63) #4:0.015(21) #5: 0.048(69) #5: 0.031(44) #5: 0.007(10) #6: 0.040(57) #6:0.025(36) #6: 0.007(10) Ave: 0.074(105) 0.070(100) Ave: 0.056(81)0.054(77) Ave: 0.042(60) 0.041(59) 3.5 cm #1: 0.080(114) #1: 0.061(87)#1: 0.047(67) #2: 0.074(106) #2: 0.060(86) #2: 0.030(43) #3: 0.061(87)#3: 0.052(74) #3: 0.009(13) #4: 0.049(70) #4: 0.039(56) #4: 0.007(10)#5: 0.041(59) #5: 0.029(41) #5: 0.007(10) Ave: 0.061(87) 0.056(80) Ave:0.048(69) 0.046(66) Ave: 0.02(29) 0.019(27)

As stated above, when the hybrid probe is compressed against the chestwall in the clinical studies, the inclusions are often located in lessthan 3 to 4 cm depth. The accuracy of reconstructed absorption maps ofthe low-contrast absorber is quite high at these depths, especially at140 Mhz. This indicates that the aforementioned method has thecapability to characterize benign lesions accurately. At a highermodulation frequency (e.g. 350 Mhz), the system detection sensitivity isreduced and the penetration of diffusive waves is also reduced.

As for the selection of the ROI or fine-mesh zone, the maximum change ofμ_(a) is 0.011 cmP^(−1P), which is obtained when the target is locatedat 2 cm and probed with 50 MHz modulation frequency (target layer #1,data not shown). This maximum change in μ_(a) will not cause any benignlesions (μ_(a)<0.15 cmP^(−1P)) to be misdiagnosed as malignant lesions(μ_(a)>0.15 cmP^(−1P)). Therefore, the region of interest selection isnot critical for reconstruction of large absorbers.

FIG. 15 shows the co-registered ultrasound images (first column ofimages) and reconstructed optical absorption maps of the 1-cmhigh-contrast absorber obtained at 50 MHz (second column of images), 140MHz (third column of images) and 350 MHz (fourth column of images)modulation frequencies, respectively. The target was translated in depthin 0.5 cm steps to centers at 1, 1.5, 2, 2.5, 3, and 3.5 cm,respectively. Table 4 provides the reconstructed average absorptioncoefficients obtained at each layer with the percentage ratio of thecalculated to the calibrated μ_(a)=0.23 cmP^(−1P) given in parenthesis.The ROI used were twice the target size (2×) and three times the targetsize (3×). The reconstructed μ_(a) has been significantly improvedacross all target depth and modulation frequencies compared to thecorresponding results for the corresponding 3 cm absorber. For 50 MHzprobing frequency, the average reconstructed values are 60%, 67%, 72%,64% 70%, and 53% for 1 cm to 3.5 cm depths, respectively, when the ROIis twice the target size. For 140 MHz, the target is reconstructed as61%, 71%, 75%, 76% 77%, and 30% for 1 cm to 3.5 cm depths, respectively,with the same ROI. The highest reconstruction accuracy is achieved at350 MHz and the reconstructed values are 74%, 79% 81% and 54%, for 1 cmto 2.5 cm target depths. No solid target mass was observed in imagesbeyond 2.5 cm depth. The better reconstruction accuracy for 350 MHz forshallow targets suggests that a higher RF modulation frequency could beuseful to probe breast lesions near the surface of the skin.

The selection of the ROI for imaging small absorbers has a significanteffect on accuracy of reconstructed target absorption coefficients.Table 4 indicates that the average reconstructed value for a 3×ROI wasapproximately 21% lower than if an ROI of twice the target size wasused. Fortunately, with the ultrasound guidance, we have the knowledgeabout the approximate target size and can select appropriate ROI fordual-mesh based optical imaging reconstruction.

TABLE 4 High-Contrast 50 Mhz 140 Mhz 350 Mhz 1-cm absorber ROI: 2 x 3 xROI: 2 x 3 x ROI 2 x 3 x 1.0 cm 0.138(60) 0.134(58) 0.141(61) 0.086(37)0.171(74) 0.099(43) 1.5 cm 0.153(67) 0.098(43) 0.163(71) 0.129(56)0.182(79) 0.129(56) 2.0 cm 0.166(72) 0.103(45) 0.173(75) 0.111(48)0.187(81) 0.131(57) 2.5 cm 0.147(64) 0.096(42) 0.175(76) 0.111(48)0.124(54) 0.083(36) 3.0 cm 0.160(70) 0.110(48) 0.177(77) 0.112(49) Notarget 3.5 cm 0.123(53) 0.098(43) 0.068(30) 0.062(27) No target

FIG. 16 shows the reconstructed optical absorption maps of the 1-cmlow-contrast absorber obtained at 50 MHz (first column), 140 MHz (secondcolumn) and 350 MHz (third column) modulation frequencies. The targetwas translated in depth in 0.5 cm steps and the target center waslocated at 1, 1.5, 2, 2.5, 3, and 3.5 cm, respectively. Becauseco-registered ultrasound images of the low-contrast absorber areessentially the same as the images of high-contrast absorber shown incolumn one of FIG. 15, these images are not shown in FIG. 16. Table 5provides the reconstructed average absorption coefficients obtained atthe target layer with the percentage ratio of the calculated to thecalibrated μ_(a)=0.07 cmP^(−1P) given in parenthesis. The low contrastsmall absorber was over-reconstructed by 0.019 cmP⁻³ on average (27%)across the depth at 50 MHz, and by 0.021 cmP⁻¹ on average (30%) at 140MHz. As discussed previously, these small changes will not cause errorsin classification of benign and malignant tumors. At 350 MHz, theover-reconstruction is 0.059 cmP^(−1P) (84%) on average. The main reasonis the higher noise level of our system at 350 MHz modulation frequency,which is more pronounced when the perturbation is smaller.

TABLE 5 Low-Contrast 50 Mhz 140 Mhz 350 Mhz 1 cm absorber ROI: 2 x 3 xROI: 2 x 3 x ROI: 2 x 3 x 1.0 cm 0.089(127) 0.059(84) 0.086(123)0.055(79) 0.114(163) 0.066(94) 1.5 cm 0.094(134) 0.063(90) 0.099(141)0.063(90) 0.139(199) 0.081(116) 2.0 cm 0.092(131) 0.063(90) 0.093(133)0.061(87) 0.134(191) 0.085(120) 2.5 cm 0.087(124) 0.065(93) 0.107(153)0.079(113) No target 3.0 cm 0.083(119) 0.032(46) 0.069(99) 0.066(94) Notarget 3.5 cm No target No targte No target

This example was conducted on a 53 year-old woman who had a 3 cm×3 cm×3cm invasive ductal carcinoma of high grade (III) located at the 6o'clock position in her right breast. The patient was scheduled forchemotherapy treatment and she was scanned her before the treatment. Twolight modulation frequencies of 50 MHz and 140 MHz were used in thescan. FIG. 17( a) shows the ultrasound image of the cancer, which islocated between 0.5 cm and 3.5 cm in depth direction. Because the cancerhas irregular margins in x-y dimensions, the ROI is chosen as the sizeof the hand-held probe (9 cm×9 cm). FIG. 17( b 1, b 2), (c 1, c 2) and(d 1, d 2) are corresponding reconstructed absorption maps obtained from50 MHz and 140 MHz at 690 nm, 780 nm and 830 nm, respectively. Similarto the 3-cm high-contrast absorber experiment, the absorption is high atthe first three target layers in the 140 MHz images. At the fourthtarget layer, the absorption mass is barely visible in 690 nm and 830 nmimages at 140 MHz. When the 50 MHz was used to probe the lesion, theabsorption mass in the fourth target layer is readily seen. Also thelight penetration at 830 nm is significantly improved at the 50 Mhzfrequency. The total hemoglobin maps and the blood oxygen saturationmaps at 50 MHz and 140 MHz are computed and the images are shown in FIG.17( e 1, e 2) and (f 1, f 2), respectively.

Table 6 lists the average absorption coefficients measured within FWHMof each target layer at three wavelengths and two modulationfrequencies, as well as computed total hemoglobin concentrations of bothmaximum and average values. The significant improvement in computedhemoglobin concentration occurs in the fourth layer and about 20(μ/liter) difference is observed. Furthermore, the estimated oxygensaturation distribution has dramatically improved because of theincreased light penetration at 830 nm. As can be seen, the oxygensaturation patterns are similar at both modulation frequencies in slice4 and 5, and the improvements occur at slice 6 and 7.

TABLE 6 50 Mhz 140 Mhz Lesion- 50 Mhz 140 Mhz 50 Mhz 140 Mhz 50 Mhz 140Mhz Hemoglobin Con. (μ/liter) Depth 690 nm 780 nm 830 nm Maximum/Ave.0.5 to #1: 0.152 0.156 #1: 0.183 0.209 #1: 0.185 0.141 125.2/86.9115.1/80.0 3.5 cm #2: 0.159 0.166 #2: 0.190 0.211 #2: 0.193 0.143125.9/91.0 112.6/81.3 #3: 0.187 0.136 #3: 0.197 0.239 #3: 0.203 0.170130.3/94.0 129.2/92.1 #4: 0.106 0.026 #4: 0.162 0.174 #4: 0.173 0.094108.7/79.9  82.1/60.7 #5: 0.043 0.026 #5: 0.085 0.064 #5: 0.106 0.014 65.2/46.1  49.3/15.3 #6: 0.027 0.026 #6: 0.044 0.019 #6: 0.073 0.014 59.6/28.2  47.7/7.7 Ave4: 0.122 0.121 Ave4: 0.183 0.208 Ave4: 0.1880.137 Ave4: 122.6/87.9 109.8/78.5 Ave5: 0.102 0.102 Ave5: 0.163 0.179Ave5: 0.172 0.111 Ave5: 111.1/79.6 597.6/65.9

It is interesting to note that the absorption maps as well as the totalhemoglobin concentration map reveal heterogeneous tumor vasculaturepatterns with more microvessels distributed at the tumor periphery. Thistype of angiogenesis distribution has been in many larger carcinomas. Inthe literature, the periphery enhancement or rim enhancement of breastlesions in contrast enhanced MR images has been reported. It istherefore concluded that the rim enhancement is due to a combination ofangiogenesis, distribution and degree of fibrosis, expression pattern ofgrowth factors, and various histologic features.

Although the higher modulation frequencies can provide improvedreconstruction accuracy, their range is limited by the reduced signalsdetected at greater penetration depths. Many components affect thesystem performance at higher RF frequencies. In the present system, theminimal bandwidth of preamplifiers, mixers, and RF switches is 500 MHz.The component that limits the frequency response is the photo multipliertube (PMT (R928)) that has 3 dB bandwidth of 160 MHz. As a result, theaverage signal decreases by approximately about 60 to about 70% from 140MHz to 350 MHz in 0.6%-0.7% Intralipid solution at a fixed PMT highvoltage gain. In addition, diffusive waves modulated at higher RFfrequencies penetrate less deep due to the exponential damping. Bothfactors affect the sensitivity of the system at the higher end of the RFmodulation (e.g. 350 Mhz). In the present system, higher coherentinterference noise at 350 MHz further limited the sensitivity andcontrast relative to the lower frequencies.

In these phantom studies, the absorbers were immersed in Intralipid™solution of 0.07% concentration. The fitted background absorption andreduced scattering coefficients were in the range of 0.02 to 0.031 cmP⁻¹and 4.9 to 6.1 cmP⁻¹ in different sets of experiments. The fitted bulkabsorption and reduced scattering coefficients from normal breasts arewithin the range of 0.01 to 0.07 cmP^(−1P) and 2 to 13 cmP⁻¹. Withhigher background tissue absorption and scattering, the penetrationdepth of diffusive waves will be further reduced.

In summary, a new frequency domain optical tomography system utilizingthree RF modulation frequencies, which are optimized for probing breastlesions of different size located at different depths was constructed.With the real-time co-registered ultrasound guidance on lesion size andlocation, an optimal light modulation frequency can be selected whichyields more accurate estimates of lesion angiogenesis and hypoxia.Phantom experiments of large and small absorbers of both high and lowoptical contrasts have demonstrated that a high modulation frequency,such as 350 MHz, is preferable for probing small lesions closer to thesurface while a low modulation frequency, such as 50 MHz, is desirablefor imaging deeper and larger lesions. A clinical example of a largeinvasive carcinoma is presented and significantmodulation-frequency-dependent and wavelength-dependent tissuepenetration of diffusive-waves has been observed.

To evaluate the possible effect of ROI or fine-mesh zone selection onaccuracy of the reconstructed absorption coefficients, a reconstructedtarget using ROIs twice the size of the real target in x-y dimensionsand three times the size of the target in x-y dimensions was used. Theresults have shown that the ROI has negligible effect on reconstructedoptical properties of larger 3-cm absorbers. For smaller 1-cm absorbers,the ROI has to be twice the size of the absorber to ensure accuratereconstruction of target optical properties. It was also found that thereadily available wavelength pair of 780 nm and 830 nm (weight(λ₁)=0.20and weight(λ₂)=0.27) is robust for estimation of total hemoglobinconcentration.

Example 4

This example was conducted to demonstrate the effects of a partiallyreflecting boundary method for diffusive optical imaging of shallowtargets in turbid media. By reflecting boundary, it is understood thatthe surface of the probe 4 in contact with tissue is non-black in color.It is desirable for the surface of the probe 4 to be reflecting i.e.,white in color.

Before focusing on imaging of shallow targets, the accuracy of usingextrapolated boundary condition to approximate partial-current boundarycondition for different reflection coefficients has to be addressed.This approximation is essential for simplifying the forward model. Thepartial-current boundary condition can be expressed formally as

$\begin{matrix}{{\varphi = {{l_{s}\frac{\partial\varphi}{\partial z}}_{z = 0}}},{{{where}\mspace{14mu} l_{s}} = {\frac{1 + R_{eff}}{1 - R_{eff}}2{D.}}}} & (10)\end{matrix}$

where φ is the fluence rate of photons and z is the direction normal tothe surface and points to the medium. R_(eff) is defined as theeffective reflection coefficient that can be determined experimentallyor by Monte Carlo simulation. D is the diffusive coefficient of photonsin tissue. According to this boundary condition, distribution of photondensity in a semi-infinite homogenous medium generated by a harmonicmodulated point source can be described by the following Green'sfunction as shown in equation (11)

$\begin{matrix}{{\varphi_{in}^{pc} = {\frac{1}{4\pi \; D}\begin{pmatrix}{\frac{\exp \left( {- {kr}_{sd}^{pc}} \right)}{r_{sd}^{pc}} + \frac{\exp \left( {- {kr}_{sid}^{pc}} \right)}{r_{sid}^{pc}} -} \\{\frac{2}{l_{s}}{\int_{0}^{\infty}\ {{l}\; {\exp \left( {{- l}/l_{s}} \right)}\frac{\exp \left( {- {kr}_{ld}} \right)}{r_{ld}}}}}\end{pmatrix}}},} & (11)\end{matrix}$

where φ_(in) ^(pc) is the fluence rate of incident wave with thepartial-current boundary condition. k is wave vector of photon densitywave. r_(sd) ^(pc) is the distance between the point source and thedetector and r_(sid) ^(pc) is the distance between the image of thepoint source and the detector. Both r_(sd) ^(pc) and r_(sid) ^(pc) arecalculated based on partial-current boundary condition. r_(ld) is thedistance between a point on the integral line l and the detector, whichis a function of the integral variable l. Equation (11) includes a lineintegral, which complicates the calculation of scattering field forimage reconstruction. Therefore, an approximated boundary conditionnamed extrapolated boundary condition has been widely used in theliterature. The corresponding Green's function is written as shown inequation (12):

$\begin{matrix}{\varphi_{in}^{ext} = {\frac{1}{4\pi \; D}\left( {\frac{\exp \left( {- {kr}_{sd}^{ext}} \right)}{r_{sd}^{ext}} - \frac{\exp \left( {- {kr}_{sid}^{ext}} \right)}{r_{sid}^{ext}}} \right)}} & (12)\end{matrix}$

where r_(sd) ^(ext) and r_(sid) ^(ext) are similar to those in Equation(11) but calculated from extrapolated boundary condition.

Based on Equations (11) and (12) and the Born approximation, thescattering fields of the photon density wave was calculated for aspherical target from two boundary conditions. The target has a radiusof 0.5 cm and is centered at 0.5 cm below the surface. The absorptioncoefficients of the target and the background are about 0.25 cm⁻¹ and0.03 cm⁻¹, respectively. The perturbation is defined as φ_(sc)^(pc)/φ_(in) ^(pc) for partial-current boundary condition and φ_(sc)^(ext)/φ_(in) ^(ext) for the extrapolated boundary condition. A relativeerror of amplitude of the perturbations between two boundary conditionsis defined as E_(am)=(|φ_(sc) ^(pc)/φ_(in) ^(pc)|−|φ_(sc) ^(ext)/φ_(in)^(ext)|)/(φ_(sc) ^(pc)/φ_(in) ^(pc)|). The phase difference ofperturbations between the two boundary conditions is defined asE_(pha)=phase(φ_(sc) ^(pc)/φ_(in) ^(pc))−phase(φ_(sc) ^(ext)/φ_(in)^(ext)). Numerical calculations have shown that if Born approximation isassumed, the E_(am) and E_(pha) are independent of the target radius andthe absorption coefficient but slightly depend on the target depth. Theerrors of amplitude and phase of perturbation E_(am) and E_(pha) as afunction of effective reflection coefficient R_(eff) are plotted inFIGS. 18 (a) and (b), respectively. The results for three typicalsource-detector separation distances (3 cm, 4 cm and 6 cm) are given inthe FIGS. 18 (a) and (b). When 0.0<R_(eff)<0.4, the amplitude error isless than 1% and phase error is less than 0.5 degrees, which arenegligible, so that the extrapolated boundary condition can beaccurately used to approximate partial-current boundary condition. WhenR_(eff)>0.7, the errors increase very quickly and become remarkable sothat the use of extrapolated boundary conditions should be avoided. When0.4<R_(eff)<0.7, the maximum errors of amplitude and phase are about 8%and 2.7 degrees, respectively, for R_(eff)=0.7 and r_(sd)=3.0 cm.

A combined probe was used with an ultrasound transducer centrallylocated to provide target size, shape, and depth information. A total of12 emitter fibers and 8 detector fibers are employed in the probe. Thedistance between emitters and detectors was distributed from about 1 cmto about 7 cm. In general, only the emitter-detector pairs withseparation distances from 2.5 cm to 6 cm can be used for imaging. Forpairs with distances less than 2.5 cm, the measured signals aresaturated; while for pairs larger than 6 cm, the measured signals aretoo weak and noisy. The entire imaging region is 9×9×3.5 (x, y and z)cm³. The dual-mesh technique discussed above was adopted for imagingreconstruction. The background region was divided into the relativelycoarse mesh and each voxel had a size of 1.5×1.5×1.0 cm³. Based on thetarget localization provided by ultrasound, the target region is limitedto a small volume of 2.5×2.5×1 cm³, which is determined by the targetsize. This inclusion region had finer voxels of 0.25×0.25×0.5 cm³ insize. In these phantom experiments, the radius of the target was fixedto about 0.5 cm and the target was located at the center of x-y plane atan approximate depth of about 0.6 to about 0.7 cm. Therefore the targetoccupied two layers near the surface. 0.6% Intralipid solution was usedto emulate soft tissue. The results obtained at 780 nm and 830 nm. Theresults obtained at 780 nm are reported below. The calibrated absorptioncoefficient of the target ranges from 0.2 to 0.3 cm⁻¹ and 0.25 cm⁻¹ wasused as an approximate value. The Born approximation was used to relatethe measured perturbation and the medium absorption distribution and aconjugate gradient technique has been used for inverse reconstruction.

The reflecting boundary was obtained by covering the black probe with analuminum foil, and its effective reflection coefficient can be estimatedby fitting the measurement data of homogeneous Intralipid to theoreticaldata of homogeneous medium obtained from the diffusion equation. It wasestimated that R_(eff)≈0.6 for a regular-aluminum foil. The amplitudeerrors were about 3.0%, 1.26% and 0.34% for emitter-detector separationdistances of 3 cm, 4 cm and 6 cm respectively. Correspondingly, thephase errors are −1.45 degrees, −1.05 degrees and −0.72 degreesrespectively.

The reconstructed phantom images are shown in the FIG. 19, in which onlythe first two target layers are shown. The slice #1 is the spatial x-yimage of the top layer, which is 0.5 cm deep in the Intralipid, and theslice #2 is the layer 0.5 cm below slice #1. FIGS. 19( a) and (b) arethe images reconstructed from an absorbing boundary (R_(eff)≈0) and areflection boundary (R_(eff)≈0.6), respectively. In FIG. 19 (a), theposition of the maximum is at (5.0, 5.0) for layer #1 and (5.0, 4.5) forlayer #2. The corresponding maximum reconstructed-values are 0.0867 and0.2967 cm⁻¹, which are 34.68% and 118.68% of the expected value (˜0.25cm⁻¹). In the FIG. 19( b), the positions of the maxima agree with theposition of the true target (4.5, 4.5). The corresponding maxima are0.212 and 0.194 cm⁻¹ for the two layers, respectively, which are 84.80%and 77.60% of the expected value. A striking difference of thereconstructed values between the two layers in the FIG. 19( a) impliesthat the absorbing boundary is not good in detecting shallow targetsless than 0.5 cm deep. However, the reflection boundary providesconsiderable uniform values in the two layers. Furthermore, in the FIG.19( a) the mean value obtained by averaging absorption coefficientswithin −6 dB of the maximum value is 0.073 for layer #1 and 0.236 cm⁻¹for layer #2. The percentages of the mean values relative to theexpected value are 29.2% for layer #1 and 94.8% for layer #2,respectively. Correspondingly, in the FIG. 19( b) they are 0.167 cm⁻¹and 0.152 cm⁻¹ and the related percentages are 66.8% and 60.8%,respectively. It is obvious that the uniformity of mean values in twolayers is improved when the reflecting boundary is used.

This work is useful in clinical applications for breast imaging.Currently, a fixed black probe is used for all patients and the shallowinclusions less than 0.5 to 0.7 cm in depth do not reconstruct correctlyin both target location and optical properties. This work provides thata probe having a reflecting surface that is not black can be used forimaging shallow inclusions

While the invention has been described with reference to exemplaryembodiments, it will be understood by those skilled in the art thatvarious changes may be made and equivalents may be substituted forelements thereof without departing from the scope of the invention. Inaddition, many modifications may be made to adapt a particular situationor material to the teachings of the invention without departing from theessential scope thereof. Therefore, it is intended that the inventionnot be limited to the particular embodiment disclosed as the best modecontemplated for carrying out this invention, but that the inventionwill include all embodiments falling within the scope of the appendedclaims.

1. A method for medical imaging of an inclusion comprising: imaging atissue volume with a probe comprising: an ultrasound transducer that isoperative to provide an on-site estimation of inclusion size andlocation; a first emitter and a first detector; the first emitter havinglight of a wavelength of about 400 to about 900 nanometers; the firstdetector detecting light of a wavelength of about 400 to about 900nanometers; a source circuit connected in operational communication tothe emitter; a detector circuit connected in operational communicationto the detector; and a central processing unit connected to the sourcecircuit and the detector circuit; scanning the tissue volume with lighthaving a wavelength of about 400 to about 900 nanometers; the tissuevolume comprising a first layer and a second layer; segmenting thescanned tissue volume into an inclusion region comprising a plurality offirst voxels and a background region comprising a plurality of secondvoxels; the volume of each second voxel being larger than the volume ofeach first voxel; identifying a tissue layer thickness and anapproximate tilting angle between the tissue layer and the probe;estimating photon density: φ(r,ω)=f(μ_(a1),μ_(s1)′,μ_(a2),μ_(s2)′),where μ_(a1), μ_(a2), μ_(s1)′, and μ_(s2)′ are absorption and reducedscattering coefficients of first and second layers, respectively;wherein absorption coefficients have the subscript “a”, scatteringcoefficients have the subscript “s”, wherein the subscript “1”represents the first layer, the subscript “2” represents the secondlayer, and the superscript (′) stands for the reduced scatteringcoefficient; applying a conjugate gradient method to estimate fittedoptical properties of the first layer and the second layer; minimizing aleast-square objective function given as min∥φ(r,ω)−φ₀(r,ω)∥²; wherein aforward Jacobian weight matrix${{Wij} = \left\lbrack {\frac{\partial\varphi_{ij}}{\partial\mu_{aj}},\frac{\partial\varphi_{ij}}{\partial D_{j}}} \right\rbrack},$that relates the photon density perturbation at detector i and imagingvoxel j with absorption coefficient change Δμ_(aj) and diffusioncoefficient change ΔD_(j), is calculated by using the bulk opticalproperties simulated from the two-layer layer model and given inequation (1) as: $\begin{matrix}{{\left\lbrack W_{ij} \right\rbrack = \begin{bmatrix}\frac{\partial\varphi_{11}}{\partial\mu_{a\; 1}} & \ldots & \frac{\partial\varphi_{1L}}{\partial\mu_{aL}} & \frac{\partial\varphi_{11}}{\partial D_{1}} & \ldots & \frac{\partial\varphi_{1\; L}}{\partial D_{L}} \\\frac{\partial\varphi_{21}}{\partial\mu_{a\; 1}} & \ldots & \frac{\partial\varphi_{2\; L}}{\partial\mu_{aL}} & \frac{\partial\varphi_{21}}{\partial D_{1}} & \ldots & \frac{\partial\varphi_{2L}}{\partial D_{L}} \\\vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\\frac{\partial\varphi_{M\; 1}}{\partial\mu_{a\; 1}} & \ldots & \frac{\partial\varphi_{ML}}{\partial\mu_{aL}} & \frac{\partial\varphi_{M\; 1}}{\partial D_{1}} & \ldots & \frac{\partial\varphi_{ML}}{\partial D_{L}}\end{bmatrix}},} & (1)\end{matrix}$ where M is the total number of detector readings; L is thetotal number of imaging voxels and φ₀(r,ω) is the photon density of thefirst layer.
 2. The method of claim 1, further comprising estimatingφ₀(r,ω), μ_(a10), μ_(a20), μ_(s10)′ or μ_(s20)′ by using a solution thathas similar values as that of the tissue.
 3. The method of claim 1,wherein the inclusion region comprising a first voxel that has a volumeof about 0.1×0.1×0.1 cubic centimeter (cm³) to about 1×1×1 cm³; thebackground region comprising a second voxel that has dimensions of about0.2×0.2×0.2 cm³ to about 2.0×2.0×2.0 cm³.
 4. The method of claim 1,wherein the scanning the tissue volume comprises scanning a cylindricalvolume of tissue having a radius large enough to approximatesemi-infinite geometry.
 5. The method of claim 1, wherein the firstlayer comprises a breast tissue, while the second layer comprises achest wall.
 6. The method of claim 1, wherein the inclusion comprises atumor, a lesion, or a cancerous region.
 7. The method of claim 1,wherein an interface between the first layer and the second layer issmooth.
 8. The method of claim 5, wherein the inclusion and a nearestsurface of the chest wall are separated by less than 2 centimeters. 9.The method of claim 1, further comprising scanning the tissue volumewith a plurality of lights, each light having a wavelength of about 400to about 900 nanometers.
 10. The method of claim 1, further comprisingscanning the tissue volume with light having a wavelength of 690, 780and 830 nanometers.
 11. The method of claim 10, wherein the light havingwavelengths of 690, 780 and 830 nanometers is emitted from laser diodes.12. The method of claim 11, wherein the laser diodes are in operativecommunication with radio frequency oscillators.
 13. The method of claim12, wherein the radio frequency oscillators operate at frequencies of350 MHz, 140 MHz and 50 MHz.
 14. The method of claim 1, wherein thescanned tissue comprises fluorescence targets or dyes or fluorophores.15. A method comprising: obtaining an image of a tissue volume with aprobe comprising: an ultrasound transducer that is operative to providean on-site estimation of inclusion size and location; a first emitterand a first detector; the first emitter having light of a wavelength ofabout 400 to about 900 nanometers; the first detector detecting light ofa wavelength of about 400 to about 900 nanometers; a source circuitconnected in operational communication to the emitter; a detectorcircuit connected in operational communication to the detector; and acentral processing unit connected to the source circuit and the detectorcircuit; scanning the tissue volume with light having a wavelength ofabout 400 to about 900 nanometers; the tissue volume comprising a firstlayer and a second layer; segmenting the scanned tissue volume into aninclusion region comprising a plurality of first voxels and a backgroundregion comprising a plurality of second voxels; the volume of eachsecond voxel being larger than the volume of each first voxel;reconstructing the image using the equation (2)[U _(sd)]_(M×1) =[W _(L) ,W _(B)]_(M×N) ×[M _(L) ,M _(B)]^(T)_(N×1)  (2) where W_(L) and W_(B) are weight matrices for the inclusionand background regions, respectively;[M_(L)] = [∫_(1_(L))Δμ_(a)(r^(′))³r^(′), …  ∫_(N_(L))Δμ_(a)(r^(′))³r^(′)]and[M_(B)] = [∫_(1_(B))Δμ_(a)(r^(′))³r^(′), …  ∫_(N_(B))Δμ_(a)(r^(′))³r^(′)]are the total absorption distributions of inclusion and backgroundregions, respectively; the weight matrices being calculated based on thebackground absorption coefficient μ _(a) and reduced scatteringcoefficient μ′_(s) measurements obtained from the normal contralateralbreast; a Born approximation being used to relate a scattered fieldU_(sd)(r_(si),r_(di),ω) measured at the source (s) and detector (d) pairi to absorption variations Δμ_(a)(r′) in each volume element of the tworegions within the sample.
 16. The method of claim 15, furthercomprising dividing the tissue volume N into N layers denoted as regionof interest #1, region of interest #2, . . . and region of interest #Nrespectively; and using the equation (2) to reconstruct one layer at atime to reconstruct the image.
 17. The method of claim 16, wherein thereconstructing each layer comprises using the Equations (3), (4) and(5):[U _(sd)]_(M×1) =[W _(L(ROI#1)) ,W _(B) ₁ ,W _(B) ₁_((other-target-layers)) ,W _(B(BG-layers)) ]×[M _(L(ROI#1)) ,M _(B) ₁,M _(B) ₁ _((other-target-layers)) ,M _(B(BG-layers)) ]T  (3)[U _(sd)]_(M×1) =[W _(L(ROI#2)) ,W _(B) ₂ ,W _(B) ₂_((other-target-layers)) ,W _(B(BG-layers)) ]×[M _(L(ROI#2)) ,M _(B) ₂,M _(B) ₂ _((other-target-layers)) ,M _(B(BG-layers)) ]T  (4). . .[U _(sd)]_(M×1) =[W _(L(ROI#N)) ,W _(B) _(N) ,W _(B) _(N)_((other-target-layers)) ,W _(B(BG-layers)) ]×[M _(L(ROI#2)) ,M _(B)_(N) ,M _(B) _(N) _((other-target-layers)) ,M _(B(BG-layers)) ]T  (5)where W_(L(ROI#1)), W_(L(ROI#2)), . . . , W_(L(ROI#N)), are weightmatrices for the region of interest #1, #2 and #N respectively, andW_(B) ₁ , W_(B) ₂ , . . . , W_(B) _(N) are weight matrices ofcorresponding background regions in the same depth layer. W_(B) ₁_((other-target-layers)), W_(B) ₂ _((other-target-layers)), W_(B) _(N)_((other-target-layers)), are weight matrices of other target layersexcluding the corresponding region of interest # target layer, andW_(B(BG-layers)) is the weight matrix of the common background layersrespectively.
 18. The method of claim 17, wherein after reconstructingeach layer, a total absorption matrix can be obtained as [M_(L(ROI#1)),M_(B) ₁ , M_(L(ROI#2)), M_(B) ₂ , . . . M_(L(ROI#N)), M_(B) _(N) ,M_(BG(BG-layers))], and an absorption distribution can be computed bydividing a corresponding voxel size in the inclusion regions and in thebackground regions.
 19. The method of claim 15, wherein the inclusionregion comprising a first voxel that has a volume of about 0.1×0.1×0.1cubic centimeter (cm³) to about 1×1×1 cm³; the background regioncomprising a second voxel that has dimensions of about 0.2×0.2×0.2 cm³to about 2.0×2.0×2.0 cm³.
 20. The method of claim 15, wherein theinclusion comprises a tumor, a lesion, or a cancerous region.
 21. Themethod of claim 15, further comprising scanning the tissue volume with aplurality of lights, each light having a wavelength of about 400 toabout 900 nanometers.
 22. The method of claim 15, further comprisingscanning the tissue volume with light having a wavelength of 690, 780and 830 nanometers.
 23. The method of claim 22, wherein the light havingwavelengths of 690, 780 and 830 nanometers is emitted from laser diodes.24. The method of claim 23, wherein the laser diodes are in operativecommunication with radio frequency oscillators.
 25. The method of claim24, wherein the radio frequency oscillators operate at frequencies of350 MHz, 140 MHz and 50 MHz.
 26. The method of claim 15, wherein thescanned tissue comprises fluorescence targets or dyes or fluorophores.27. The method of claim 15, wherein the first emitter or the firstdetector is inclined at an angle of about 1 to about 30 degrees withrespect to a surface of the probe that contacts tissue.
 28. The methodof claim 15, further comprising improving the accuracy of thereconstructed image using the Equation (6)[U _(sd)]_(M×1) =[W _(L) ,W _(B) ,W _(S)]_(M×N) ×[M _(L) ,M _(B) ,M_(S)]^(T) _(N×1)  (6) where W_(L) and W_(B) are absorption weightmatrices for lesion and background regions, respectively; W_(S), is thescattering weight matrix for the entire imaging volume;[M_(L)] = [∫_(1_(L))Δμ_(a)(r^(′))³r^(′), …  ∫_(N_(L))Δμ_(a)(r^(′))³r^(′)]   and  [M_(B)] = [∫_(1_(B))Δμ_(a)(r^(′))³r^(′), …  ∫_(N_(B))Δμ_(a)(r^(′))³r^(′)]are the total absorption distributions of inclusion and backgroundregions, respectively.[M_(S)] = [∫₁Δ D(r^(′))³r^(′), …  ∫_(N)Δ D(r^(′))³r^(′)]is the total diffusion distribution of the entire tissue volume; andwherein the respective weight matrices are calculated based on thebackground absorption coefficient μ _(a) and reduced scatteringcoefficient μ′_(s) measurements obtained from a normal contralateralbreast.
 29. The method of claim 1, further comprising estimatingdeoxygenated (deoxyHb), oxygenated (oxyHb) hemoglobin and totalhemoglobin concentration totalHb(r′)=deoxyHb(r′)+oxyHb(r′) andoxygenation saturation${Y\; \%} = {\frac{{{oxy}{Hb}}\left( r^{\prime} \right)}{{{{oxy}{Hb}}\left( r^{\prime} \right)} + {{deoxy}\left( r^{\prime} \right)}}100\%}$from Equations (7), (8) and (9) respectively: $\begin{matrix}{{{{total}{Hb}}\left( r^{\prime} \right)} = {{\frac{1}{\Delta}\left\{ {{\left\{ {ɛ_{{HbO}_{2}}^{\lambda_{2}} - ɛ_{Hb}^{\lambda_{2}}} \right\} {\mu_{a}^{\lambda_{1}}\left( r^{\prime} \right)}} + {\left\{ {ɛ_{Hb}^{\lambda_{1}} - ɛ_{{HbO}_{2}}^{\lambda_{1}}} \right\} {\mu_{a}^{\lambda_{2}}\left( r^{\prime} \right)}}} \right\}}\mspace{146mu} = {{{{weight}\left( \lambda_{1} \right)}{\mu_{a}^{\lambda_{1}}\left( r^{\prime} \right)}} + {{{weight}\left( \lambda_{2} \right)}{\mu_{a}^{\lambda_{2}}\left( r^{\prime} \right)}}}}} & \begin{matrix}(7) \\(8)\end{matrix} \\{{{{and}\mspace{14mu} Y\%} = {\frac{{{- ɛ_{Hb}^{\lambda_{2}}}\frac{\mu_{a}^{\lambda_{1}}\left( r^{\prime} \right)}{\mu_{a}^{\lambda_{2}}\left( r^{\prime} \right)}} + ɛ_{Hb}^{\lambda_{1}}}{{\left( {ɛ_{{HbO}_{2}}^{\lambda_{2}} - ɛ_{Hb}^{\lambda_{2}}} \right)\frac{\mu_{a}^{\lambda_{1}}\left( r^{\prime} \right)}{\mu_{a}^{\lambda_{2}}\left( r^{\prime} \right)}} - \left( {ɛ_{{HbO}_{2}}^{\lambda_{1}} - ɛ_{Hb}^{\lambda_{1}}} \right)}100\%}},} & (9)\end{matrix}$ where μ_(a) ^(λ) ¹ (r′) and μ_(a) ^(λ) ² (r′) are opticalabsorption at volume element r′ and Δ=ε_(Hb) ^(λ) ¹ ε_(HbO) ₂ ^(λ) ²−ε_(HbO) ₂ ^(λ) ¹ ε_(Hb) ^(λ) ² .
 30. An apparatus for medical imagingusing diffusive optical tomography and fluorescent diffusive opticaltomography comprising; a probe comprising a first emitter and a firstdetector; the first emitter or the first detector being inclined at anangle θ of up to about 30 degrees to a surface of the probe thatcontacts tissue; the probe comprising a an ultrasound transducer; asource circuit connected in operational communication to the emitter; adetector circuit connected in operational communication to the detector;a central processing unit connected to the source circuit and thedetector circuit; a display operably connected to the central processingunit; and, wherein the central processing unit is capable of processinginformation to provide diffusive optical tomography and fluorescentdiffusive optical tomography.
 31. The apparatus of claim 30, wherein theultrasound transducer emits energy having a frequency of 20,000 to abillion hertz.
 32. The apparatus of claim 30, wherein the probe has areflecting surface that contacts tissue.
 33. The apparatus of claim 30,wherein the probe has a white surface.
 34. The apparatus of claim 30,wherein the probe has a black surface that absorbs substantially all ofthe light incident upon it.